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Abstract 


This  Phase  I  Strategic  Technology  Transfer  Report  (STTR)  final  report  summarizes  the 
development  of  a  prototype  interactive  software  environment,  SpecLab,  for  analyzing  non¬ 
stationary  processes  that  admit  local  power-law  representations.  Standard  spectral  analysis 
procedures  assume  stationarity,  whereas  most  naturally  occuring  processes  admit  random 
departures  from  strict  stationarity.  SpecLab  estimates  and  synthesizes  the  non-stationary 
process  by  allowing  both  the  power-law  parameters  and  the  power-law  scale  range  to  vary 
over  a  data  segmentation  chosen  interactively  by  the  user.  The  SpecLab  procedures  are 
accessible  via  a  graphical  user  interface  that  guides  the  user  through  the  steps  involved  in 
selecting  data  segmentations  and  executing  the  estimation  procedures. 

In  its  final  form  SpecLab  will  provides  efficient  user  access  to  leading-edge  analysis 
procedures  for  non-stationary  processes.  SpecLab  is  also  configured  to  provide  reproducible 
research  that  would  ordinarily  be  available  only  as  text,  equations,  and  graphs.  The  Phase 
1  prototype  demonsrates  all  the  essential  features  of  a  marketable  software  product  to  be 
developed  under  a  Phase  II  effort. 

The  report  is  organized  as  follows.  Chapter  1  provides  introductory  background  infor¬ 
mation  and  summarizes  the  developmental  Phase  I  accomplishments  against  the  content 
of  the  original  plan.  Chapter  2  provides  a  detailed  functional  description  of  the  prototype 
software  that  comprises  SpecLab  Version  1.  Chapter  3  summarizes  the  proposed  Phase  II 
STTR  effort  that  will  develop  the  final  SpecLab  product.  Chapter  4  reproduces  the  paper 
“Wavelet  Based  Estimation  of  Local  Kolmogorov  Turbulence,”  by  G.  Papanicolaou  and  K. 
S0lna,  which  has  bee  submitted  for  publication.  The  paper  summarizes  the  theoretical  un¬ 
derpinnings  of  SpecLab  and  demonstrates  the  reproducible  research  capability.  Chapter  5 
reviews  the  relation  of  standard  spectral  analysis  procedures  and  wavelet-based  estimation. 
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Chapter  1 

SpecLab  Overview 


1.1  Introduction 

Stochastic  models  that  incorporate  non-stationary  structure  through  long-range  correlations 
have  been  the  focus  of  much  attention.  Such  models  have  been  used  to  analyze  communi¬ 
cations  data  involving  local  area  networks,  video,  wide  area  networks,  and  common  channel 
signaling.  In  geophysics  long-range  dependencies  have  been  found  in  the  fluctuations  of 
parameters  that  characterize  the  earth’s  crust.  Finally,  stochastic  models  that  incorporate 
long-range  dependence  are  fundamental  in  the  analysis  of  turbulent  flow  data.  This  latter 
application  provided  the  focus  of  and  motivation  for  developing  SpecLab. 

The  theory  of  stationary  stochastic  processes  is  well  developed,  but  deviations  from  sta¬ 
tionary  are  crucial  in  the  aforementioned  applications  and  in  many  other,  if  not  most,  ap¬ 
plications  that  involve  stochastic  processes.  Yet,  modeling  and  estimation  of  non-stationary 
processes  is  neither  well  developed  nor  an  integral  part  of  mainstream  data  analysis  proce¬ 
dures.  SpecLab  makes  extensive  use  of  the  wavelet  transform,  to  provide  spatially  localized 
information  about  the  essential  frequency  content  of  the  data.  An  important  objective  of  the 
SpecLab  development  is  to  make  this  recent  progress  in  the  development  of  practical  meth¬ 
ods  for  dealing  with  non-stationary  long-range  dependent  processes  available  and  accessible 
to  a  broad  spectrum  of  potential  users. 

To  achieve  this  objective,  we  have  structured  SpecLab  as  an  interactive  software  envi¬ 
ronment  for  data  analysis,  estimation  and  synthesis  of  local  power-law  processes,  which  are 
described  in  more  detail  below  and  in  Chapter  4.  Local  power-laws  provide  a  robust,  well 
understood  model  for  the  development.  The  interactive  structure  of  SpecLab  allows  data 
exploration  by  identifying  a  representative  power-law  structure  and  synthesizing  realizations 
of  the  data  from  the  estimated  model  parameters.  SpecLab  also  serves  as  a  learning  en¬ 
vironment  for  fractals,  multifractals  and  local  power  laws  in  general.  Finally,  it  serves  as 
a  tool  for  reproducible  research  in  the  subject  area.  Reproducible  research  is  intended  to 
make  the  data  and  the  analysis  procedures  that  underlie  research  reports  available  along 
with  the  textual  and  graphic  descriptors. 

1.2  Local  power-law  processes 

Fractional  Brownian  motion  (fBm)  represents  a  local  power-law  process  in  its  simplest  form. 
Fractional  Brownian  motion,  Bn ,  is  a  non-stationary  Gaussian  process  with  stationary  in¬ 
crements.  The  variance  of  the  stationary  increment  process  is  characterized  by  the  structure 
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function 


E[(B„{x)  -  BH(x  -  Ax))2]  =  a2\Ax\2H.  (1.1) 

The  Hurst  exponent  H  determines  the  correlation  distance  for  the  increments  of  the  pro¬ 
cess  while  cr 2  establishes  the  absolute  level  of  the  correlation.  Ordinary  Brownian  motion 
corresponds  to  H  =  1/2  and  its  consecutive  increments  can  be  shown  to  be  independent. 
Increasing  and  decreasing  values  of  the  Hurst  exponent  correspond  to  positive  and  negative 
correlations,  respectively  (see  Section  4.2).  Fractional  Brownian  motion  is  self-similar  since 
Bh{x )  and  aH Bn(x/a)  have  the  same  distributions.  The  process  is  essentially  invariant 
with  respect  to  scaling  and,  thus,  contains  structural  features  on  all  scales.  This  attribute 
reflects  the  long-range  dependence  aspect  of  the  process.  Although  the  structure  function 
(1.1)  has  a  power-law  form,  the  complementary  spectral-domain  power  law  is 

the  source  of  the  definition.  The  spectral  strength  C  is  related  to  the  Hurst  exponent  and 
a  (see  Chapter  5). 

The  local  power-law  model  used  in  SpecLab  generalizes  the  fBm  model  in  two  ways. 
First,  the  power  law  parameters  a  and  H  are  not  constant,  but  vary  slowly  with  respect 
to  location  x.  These  variations  are  modeled  as  secondary  stochastic  processes.  Second,  the 
model  (1.1)  is  applied  only  over  a  subset  of  scales:  Xin  <  x  <  xout,  which  is  called  the 
inertial  range  in  turbulence  theory.  Turbulence  theory  is  an  important  application  of  the 
fBm  model.  Kolmogorov’s  scaling  law  yields  a  spectral  index  for  temperature  fluctuations 
that  corresponds  to  H  =  1/3  in  the  inertial  range.  However,  real  data  typically  show  a 
multifractal  behavior  as  manifest  by  variations  in  a  and  H  as  well  as  variations  in  the 
inertial  range  itself.  SpecLab  is  designed  to  estimate  these  variations. 

Selection  of  an  appropriate  bias- variance  trade-off  is  critical  to  the  estimation  procedure. 
For  example,  if  the  mean  power-law  parameters  are  estimated  over  segments  of  fixed  length, 
the  fluctuations  in  the  estimate  decrease  with  the  length  of  the  segments;  however,  the 
bias  increases  with  the  length  of  the  segments  since  intrinsic  variations  in  the  power-law 
parameters  are  smoothed  out.  The  local  power-law  model  and  the  scheme  designed  to 
handle  the  above  estimation  problem  are  described  in  more  detail  in  Chapter  4. 


1.3  Summary  of  Phase  I  STTR  Effort 

In  our  Phase  I  STTR  proposal  the  following  classes  of  analysis  routines  were  identified  as 
the  functional  units  of  a  viable  software  environment  for  analyzing  non-stationary  power-law 
processes. 

Data  analysis  and  simulation 

These  are  basic  routines  that  transform  the  data,  perform  wavelet  scale-filtering,  and 
display  the  data  in  various  ways.  This  set  of  procedures  also  contains  routines  that 
simulate  local  power-law  processes.  The  most  elaborate  of  these  routines  is  the  fBm 
simulator. 

Scale  spectra  estimation 

These  are  routines  that  analyze  the  transformed  data  obtained  from  the  basic  analysis 
and  simulation  procedures.  The  collection  includes  routines  that  estimate  scale  spec¬ 
tra,  the  inertial  range,  and  the  power-law  parameters.  Scale  spectra  are  described  in 
4.2.2.  Note  that  the  power-law  parameter  estimates  conain  an  estimation  error  that 
depends  on  the  segmentation  of  the  data  vector.  This  error  is  approximately  removed 
by  the  filter  described  in  Section  4.4. 
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Segmentation  analysis 

The  first  group  of  routines  in  this  set  is  designed  to  extract  an  optimal  smoothing 
of  the  power  law  parameters  derived  from  the  scale  spectra  associated  with  a  given 
segmentation.  The  smoothing  is  based  on  an  estimate  of  the  spatial  correlation  of  these 
parameters  over  segments,  as  described  in  Section  4.4.3.  SpecLab  includes  tools  that 
support  the  estimation  of  the  structure  function  that  defines  the  spatial  correlation 
for  a  given  power-law  process  realization. 

The  second  group  of  routines  in  this  set  provides  guidelines  for  choosing  an  optimal 
prior  segmentation  of  the  signal.  The  optimal  segmentation  balances  bias  and  variance 
as  discussed  in  the  introductory  description. 

Fourier  spectral  estimation 

For  completeness  and  as  a  means  of  relating  the  SpecLab  analysis  to  results  obtained 
from  conventional  analyses,  standard  spectral  analysis  routines  are  included.  The 
relations  are  summarized  in  Chapter  5. 

Interactive  interface  and  intermediate  level  routines 

This  is  a  graphical  user  interface  (GUI)  that  allows  the  user  to  apply  SpecLab  routines 
in  a  user-friendly  way.  Intermediate-level  routines  that  carry  out  complete  subtasks 
like  optimal  smoothing  of  spectral  parameters  are  invoked  by  the  GUI.  The  interme¬ 
diate  level  routines  make  it  convenient  for  the  user  to  identify  the  key  computational 
procedures  in  SpecLab  and  how  they  can  be  tailored  toward  specialized  applications. 

Demos 

These  are  Matlab  demonstration  scripts  that  perform  end-to-end  analyses  of  specific 
data  sets.  They  serve  as  an  introduction  to  SpecLab  and  provide  the  SpecLab  user 
with  practice  exercises  to  explore  its  capabilities. 

Our  Phase  I  plan  has  been  successfully  completed.  The  SpecLab  prototype  software 
provides  an  interactive  Matlab-based  software  environment  for  analysis  of  local  power-law 
processes.  Users  can  select  sample  data  sets,  generate  synthetic  data  sets,  or  import  data 
of  their  own  choosing.  SpecLab  procedures  will  display  the  input  data,  its  spectrum,  and 
various  parameters  that  will  assist  in  the  subsequent  processing.  The  package  is  designed  to 
provide  the  user  with  a  systematic  method  for  identifying  appropriate  inertial  ranges  and  the 
intrinsic  variation  in  the  power-law  parameters.  The  SpecLab  prototype  developed  under  the 
Phase  I  STTR  uses  Matlab’s  GUI  facilities  and  Matlab ’s  system  for  on-line  documentation 
and  help.  It  requires  both  Matlab  and  the  Wavelet  Tool  Box,  but,  under  the  proposed  Phase 
II  effort,  stand-alone  Unix  and  PC  versions  will  be  developed. 

The  SpecLab  estimation  procedures  have  been  demonstrated  to  work  well  both  with  real 
and  synthetic  data.  Moreover,  all  of  the  results  discussed  in  Section  4  can  be  reproduced 
with  SpecLab,  which  demonstrates  its  reproducible  research  capability.  SpecLab  has  been 
implemented  with  a  GUI.  This  task  was  deemed  sufficiently  important  that  we  accelerated 
the  development  beyond  the  original  Phase  I  plan. 

Functions  that  have  not  yet  been  implemented  include  (i)  Adaptation,  which  is  a  pro¬ 
cedure  to  help  the  user  select  a  prior  segmentation,  (ii)  Measurement  noise,  which  is  a 
procedure  that  quantifies  measurement  noise  in  the  data,  and  (iii)  Spectral  comparison.  The 
material  described  in  Chapter  5  covers  (iii).  The  incorporation  of  these  remaining  functions 
is  discussed  in  Chapter  3. 
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Chapter  2 


SpecLab  Version  I 


2.1  Introduction 

SpecLab  is  a  Matlab-based  prototype  software  environment  designed  to  illustrate  the  es¬ 
sential  elements  of  an  interactive  analysis  tool  for  non-st ationary  processes.  It  is  based  on 
the  work  done  by  G.  Papanicolaou  and  K.  S0lna  with  the  collaboration  of  V.  Kruger  and 
C.  Rino  of  Vista  Research  and  D.  Washburn  of  the  Air  Force  Lab  at  Kirtland  Air  Force 
Base.  SpecLab  was  developed  and  written  by  K.  Solna.  The  software  as  described  below 
is  available  via  anonymous  ftp  upon  written  request  to  the  authors.  The  remainder  of  this 
section  takes  the  form  of  a  preliminary  user’s  guide  for  SpecLab  Phase  I. 


2.2  Running  SpecLab 

SpecLab  is  written  in  interactive  Matlab  version  5.3.  This  user’s  guide  assumes  that  the  user 
is  running  Matlab  and  is  familiar  with  its  operation  and  data  structures.  It  is  also  necessary 
to  have  the  Matlab  Wavelet  Toolbox  installed.  SpecLab  will  run  under  the  Unix/Linux  or 
the  PC/NT  versions  of  Matlab.  The  startup  script  will  define  the  appropriate  environment. 
In  the  remainder  of  this  section  we  assume  that  the  SpecLab  directory  has  been  downloaded 
and  uncompressed. 

With  Matlab  running,  change  directories  to  ‘SpecLab/Lab’  to  set  up  the  program  en¬ 
vironment.  The  procedure  is  the  same  for  Linux  or  PC  versions.  To  start  SpecLab  type 
‘SpecLab’. 
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2.2.1  Reading  or  simulating  a  data  vector 

After  typing  SpecLab  in  Matlab  the  user  will  be  presented  with  the  following  menu: 


CHOOSE  INPUT: 


HELP 


BINARY  file  with  data- vector 


ASCII  file  with  data  vector 


SIMULATE  a  process 


DEMO 


QUIT 


The  first  option  initiates  an  on-line  help  utility.  The  second  and  third  options  load  existing 
binary  or  ASCII  files.  The  fourth  option  allows  the  user  to  synthesize  a  data  vector  with 
prescribed  local  power-law  characteristics.  The  demo  option  will  walk  the  user  through  the 
various  processing  steps  to  illustrate  the  functional  elements  of  SpecLab. 

For  the  second  or  third  options  a  menu  appears  that  identifies  the  following  precomputed 
data  vectors: 

•  pure  Brownian  motion 

•  Brownian  motion  with  tapered  spectrum 

•  a  pure  Kolmogorov  process  with  H  =  1/3 

•  local  Brownian  motion 

•  typical  aerothermal  data  (the  data  set  described  in  Chapter  4) 

The  user  can  alternatively  specify  the  name  of  a  file  that  contains  data  stored  as  a  sequence 
of  numbers  in  a  one-column  format.  The  routines  implementing  file  communication  are 
located  in  subdirectory  ‘SpecLab/fileio.’ 

The  simulation  option  initiates  a  sequence  of  menus  that  input  the  following  local  power- 
law  parameters: 

•  length  of  data  vector 

•  inner  scale 

•  outer  scale 

•  correlation  length  of  parameter  variations 

•  mean  Hurst  exponent 

•  fluctuations  of  Hurst  exponent  H 

•  fluctuations  of  the  log  intercept  log(cr2) 

The  parameters  H  and  a2  can  be  chosen  to  be  deterministic  or  random. 

Finally,  choosing  the  demo  option  initiates  a  workout  that  illustrates  the  estimation  of 
a  local  power-law  process,  as  described  in  Section  2.5. 
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2.2.2  Analyzing  and  plotting  data 

Once  SpecLab  has  loaded  or  synthesized  a  data  vector,  the  user  is  presented  with  a  menu  that 
identifies  the  various  procedures  that  can  be  performed  on  the  data  vector  as  summarized 
below: 


CHOOSE  INPUT: 

HELP 

CHOOSE  segmentation 

ESTIMATE  inertial  range 

ESTIMATE  model  for  power  law  parameters 

PLOT  process 

PLOT  wavelets 

PLOT  scale  spectrum  (color) 

PLOT  scale  spectrum  (graphs) 

PLOT  power  law  parameters 

PLOT  smoothed  power  law  parameters 

PRINT  current  plot 

SAVE  current  data 

FILTER  data 

RETURN  to  input  menu 

QUIT 
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As  a  guide  to  the  choices,  recall  that  estimating  the  local  power-law  model  parameters 
involves  three  steps: 

(i)  Segmentation  of  the  raw  data  vector  into  a  specified  number  of  equal-length  segments. 

(ii)  Estimation  of  inertial  range  for  each  segment  and  of  the  power-law  parameters  for  each 

segment. 

(iii)  Estimation  of  the  structure-function  model  that  characterizes  the  correlation  of  the  of 
the  power-law  parameters  for  the  segmentation. 

Selecting  the  second,  third,  and  fourth  entries  in  sequence  carries  out  these  operations. 
The  main  menu  also  has  a  number  of  options  for  plotting  the  data  that  can  be  executed 
before  and  during  the  power-law  parameter  estimation.  The  action  choices  are  described  in 
more  detail  below. 

(i)  Segment  the  data 

Selecting  the  segmentation  option  initiates  a  secondary  menu  with  the  following  choices: 

(a)  Enter  the  number  of  equal-length  segments. 

(b)  Display  the  space-scale  spectra  in  a  color  format. 

(c)  Display  the  scale  spectra  as  graphs. 

(d)  Return  to  main  menu  with  current  segmentation. 

The  scale  spectra  are  calculated  over  each  segment  and  stored  in  the  SpecLab  data 
structure;  see  Section  2.3.  Once  this  step  has  been  performed,  the  scale-spectra  displays 
can  be  initiated  from  the  main  menu.  The  default  segmentation  is  one  segment. 

(ii)  Estimate  inertial  range 

Selecting  this  estimation  option  initiates  a  secondary  menu  with  the  following  choices: 

(a)  Estimate  the  inertial  range. 

(b)  Specify  the  inertial  range. 

(c)  Display  the  space-scale  spectra  in  a  color  format. 

(d)  Display  the  scale  spectra  as  graphs. 

(e)  Return  to  main  menu  with  current  segmentation. 

For  option  (a)  the  interial  range  estimates  are  computed  automatically  based  on  com¬ 
parison  of  the  scale  spectra  with  pure  fractional  Brownian  motion;  see  Section  4.4.2.  For 
option  ( b )  the  user  chooses  a  constant  segmentation  that  is  applied  to  all  data  segments. 
The  default  inertial  range  uses  all  available  scales. 

(iii)  Estimate  structure  function 

The  stochastic  variation  of  the  power-law  parameters  from  segment  to  segment  is 
modeled  by  independently  fitting  an  exponential  structure  functions  of  the  form 

w  +  v(l  —  exp(— A  x/l)) 

to  H  and  log(a2).  Here  l  is  the  ‘correlation  ranged  w  is  the  ‘intercept,’  v  is  the  ‘vari¬ 
ance’  and  Ax  converts  the  lag  number  to  sample  units.  The  user  is  supported  in 
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choosing  the  appropriate  structure  function  with  displays  that  compare  the  empiri¬ 
cal  structure  function  calculated  from  the  data  with  the  default  model  (exponential) 
structure  function.  The  user  can  modify  the  parameters  in  the  model  interactively  in 
order  to  improve  the  fit  via  a  secondary  menu  with  the  following  options: 

•  Change  correlation  range. 

•  Change  variance. 

•  Change  intercept. 

•  Change  max  variogram  lag. 

•  Accept  current  parameters. 

Note  that  if  the  number  of  segments,  and  hence  the  length  of  the  power-law  parameter 
processes,  is  too  short,  then  the  default  structure  function  model  is  used.  The  default 
model  corresponds  to  constant  power-law  parameters,  which  is  also  assumed  if  the 
estimation  of  structure  function  step  is  omitted. 

(iv)  Scale  filter  data 

Scale  filtering  is  a  process  by  which  the  data  are  reconstructed  using  only  a  selected 
range  of  the  scale  spectra  values.  The  operation  is  useful  as  a  preprocessing  operation 
to  remove  spurious  contamination  of  real  data  vectors. 

Choosing  this  option  brings  up  a  secondary  menu  with  the  following  options: 

•  Select  scale  range. 

•  Display  scale-filtered  process. 

•  Return  with  filtered  data  replacing  original  data. 


2.3  SpecLab  data  structure 

The  data  vector  processed  by  SpecLab  is  a  data  structure  that  accumulates  inputs  as  the 
various  operations  are  performed.  The  scripts  in  SpecLab  communicate  using  this  data 
structure.  An  example  of  a  data  structure  is  given  below 

Z  = 


data: 
wavec : 
wavel : 
nseg : 
scspect : 
inertial: 
parproc : 
f iltpar : 
f iltparproc : 
trueparproc : 


[65536x1  double] 

[65536x1  double] 

[12x1  double] 

32 

[32x10  double] 

[32x2  double] 

[32x2  double] 

[0.0026  0  7.  0.0463  08.] 
[32x2  double] 

[] 


The  field  values  (see  also  Section  5.2.2)  contain  the  following  objects: 
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data 

wavec 

wavel 

nseg 

scspect 

inertial 

parproc 

filtpar 

filtparproc 

trueparproc 


the  data  vector 

the  vector  c  in  Matlab’s  wavelet  decomposition 

the  vector  1  in  Matlab’s  wavelet  decomposition 

number  of  segments  in  the  uniform  segmentation  of  the  data  vector 

the  scale  spectra  for  the  given  data  vector  and  segmentation 

the  estimated  inertial  ranges  for  the  given  segmentation 

the  power-law  parameter  estimate  based  on  the  scale  spectra 

and  estimated  inertial  ranges 

the  parameters  in  the  structure-function  models  for  the  power 
law  parameters 

the  filtered  power-law  parameters 

if  the  data  vector  was  simulated  the  true  parameter  processes  are 
stored  here 


2.4  SpecLab  directory  structure 

The  SpecLab  directories  and  files  are 

Lab  InitZ  Main  SpecLab 

demoAero  demoBm  demolnert  demoKol  demoLfBm 
specpath  startup 

display  displayP  displays  displaySyn  display Vario  displayW  displayX 

fileio  ReadSpec  WriteSpec  getdat  writedat 

SpecLab 

options  Helplnput  HelpSpec  InertialSpec  ProcSpec  ScaleSpec 
ScfiltSpec  SmoothSpec  SynthSpec  VarioSpec 

demo  aero  pause  showinput  showmenu 

tools 

inertialest 

Getlnertial  Inertial  MakeRefres  Refres 
vario  est 

Get  Vario  VarioEst  regmodel 
spectest 

GetParproc  Linreg  Spect  hfunc 

filter 

Filt_Par  Scale_Filt  Taper 

simulation 

SegSim  SimExp  sim_ob  simcoef  vfunc 
The  scripts  are  documented  in  Appendix  A. 
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2.5  The  “Workout” 


The  workout  shows  how  SpecLab  commands  can  be  used  to  generate  all  the  figures  in 
Section  4.  The  script  that  implements  the  workout  is  ‘aero.m’  and  is  located  in  directory 
‘SpecLab/aero.’  The  demo  illustrates  the  steps  involved  in  estimating  a  local  power-law 
model.  By  looking  at  the  script  itself  the  user  can  get  information  about  how  the  interme¬ 
diate  level  SpecLab  routines  implement  the  estimation  steps.  The  user  starts  the  workout 
by  choosing  the  option  ‘demo’  in  the  input  menu  described  in  Section  2.2.1.  Starting  the 
demo  gives  the  following  figures 

Figure  4.1  Displays  the  data  vector.  The  data  vector  is  ‘demoAreo’  located  in  directory 
‘SpecLab/Lab.’  This  data  vector  can  be  analyzed  independently  by  choosing  fhte  second  or 
third  options  in  the  input  menu  described  in  Section  2.2.1. 

Figure  4.2.  Next,  the  user  is  presented  with  the  main  menu  given  in  Section  2.2.2. 
Choosing  the  menu  option  ‘PLOT  wavelets’  creates  the  figure. 

Figure  4.3.  The  user  is  asked  to  choose  the  main  menu  option  ‘CHOOSE  segmentation’ 
and  then  give  the  number  of  segments.  Choosing  the  ‘plot  scale  spectra  (graphs)’  option 
creates  the  figure. 

Figure  4.7.  The  user  is  asked  to  choose  the  main  menu  option  ‘ESTIMATE  inertial 
range’  and  then  choose  the  ‘estimate  inertial  range’  option.  Choosing  the  ‘plot  inertial  scale 
spectra  (graphs)’  option  creates  the  figure. 

Figure  4.8.  By  choosing  the  main  menu  option  ‘PLOT  Power  law  parameters’  the 
user  obtains  the  figure  with  the  estimated  parameter  processes.  Note  that  the  parameter 
processes  have  been  truncated  so  that  only  the  parameters  in  the  latter,  high-turbulence 
regime  are  given. 

Figure  4.10.  Finally,  by  choosing  the  main  menu  option  ‘PLOT  Smoothed  power  law 
parameters’  the  user  obtains  the  figure  with  smoothed  parameter  processes.  Before  doing 
the  smoothing  of  the  parameter  processes  the  structure  function  for  the  parameter  processes 
need  to  be  estimated.  In  the  workout  these  were  chosen  to  be  the  ones  used  in  [24],  when 
analyzing  another  data  set  these  must  be  estimated  choosing  the  ‘ESTIMATE  model  for 
power-law  parameters’  option  in  the  main  menu. 
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Chapter  3 


Phase  II  Plan 


3.1  Introduction 

The  prototype  SpecLab  environment  demonstrates  all  the  essential  elements  of  an  interactive 
analysis  tool  for  non-stationary  local  power-law  processes.  In  its  present  form,  however,  it 
has  a  number  of  limitations  that  would  preclude  its  use  by  all  but  knowledgeable  users. 
Under  a  Phase  II  effort  SpecLab  would  be  converted  into  a  more  general  and  complete 
software  environment,  together  with  Web-based  documentation  and  interactive  capability. 
The  general  structure  of  SpecLab  lends  itself  readily  to  these  enhancements. 

3.1.1  Analysis  and  synthesis  of  local  power  laws 

We  have  emphasized  that  while  the  theory  of  stationary  stochastic  processes  is  well  devel¬ 
oped,  there  are  few  widely  available  utilities  for  analyzing  non-stationary  process.  Non- 
stationary  processes  by  their  very  nature  require  adaptive  processing.  Thus,  we  believe  that 
one  of  the  most  important  features  of  SpecLab  is  its  capability  for  interactive  and  adaptive 
processing.  This,  in  turn,  places  more  critical  demands  on  user  support.  SpecLab  ultimately 
will  make  these  techniques  available  in  a  robust  and  user-friendly  software  environment.  In¬ 
teractive  analysis  and  synthesis  for  identification  and  replication  of  non-stationary  effects 
in  a  large  class  of  fractal  and  multifractal  processes  are  novel  features  that  SpecLab  will 
provide.  SpecLab  also  contains  a  number  of  features  that  extend  beyond  those  associated 
with  local  power  laws.  With  the  proposed  software  environment  these  features  can  be  used 
for  more  general  problems.  Thus,  the  focus  on  local  power-laws  is  not  overly  restrictive. 

3.1.2  Reproducible  Research 

Reproducible  research  was  pioneered  in  the  development  of  the  software  package  Wavelab  [4]. 
The  stated  objective  was  to  make  the  full  content  of  research  available  via  a  complete  soft¬ 
ware  environment  that  implements  the  computational  procedures  used  to  generate  reported 
numerical  results.  The  development  of  SpecLab  will  fully  preserve  this  capability.  While 
the  SpecLab  reproducible  research  umbrella  is  specific  to  the  analysis  of  non-stationary  pro¬ 
cesses  using  a  local  fBm-based  power-law  model,  the  structure  is  a  template  for  distinctly 
different  applications  that  are  based  on  interpretive  analysis  tools. 
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3.1.3  SpecLab  as  a  learning  tool 

Fractal  and  multifractal  processes  have  been  the  subject  of  much  recent  interest.  SpecLab 
will  serve  as  a  learning  tool  for  users  who  want  to  understand  such  processes.  Through  that 
understanding  we  expect  new'  applications  to  emerge.  A  key  aspect  of  SpecLab  is  analysis 
and  processing  using  the  wavelet  transform.  The  package  contains  user  friendly  procedures 
for  filtering  in  scale  and  time  and  visualization  of  the  considered  process  based  on  this 
transform.  Thus,  SpecLab  can  also  serve  as  an  introduction  to  wavelet  applications. 

3.1.4  Creating  a  portable  software  environment 

The  creation  of  the  world-wide  Web  has  revolutionized  information  exchange.  We  will  make 
use  of  the  Web  to  provide  easy  access  to  SpecLab  and  its  documentation.  By  making  the 
product  easily  available  in  a  platform-independent  form,  the  impact  of  research  and  devel¬ 
oped  software  is  enhanced  enormously.  We  expect  that  the  development  of  this  aspect  of 
the  software  package  will  formalize  an  emerging  trend  at  universities  and  research  institutes. 

3.2  Project  outline 

Under  the  Phase  II  effort  the  utilities  described  below  would  be  developed. 

3.2.1  Extended  user  guide 

•  SpecLab  Web  page.  A  Web  page  will  be  developed  that  will  give  background 
information  about  SpecLab  and  its  functionality.  It  will  be  based  on  an  updated 
version  of  the  information  in  this  Phase  I  report  for  SpecLab.  That  is,  the  Web 
page  will  identify  the  SpecLab  objectives,  give  information  about  locally  stationary 
processes,  and  document  the  software  along  the  lines  of  an  enhanced  HTML  version  of 
Chapter  2.  The  visitor  to  the  Web  page  will  also  be  able  to  download  relevant  papers, 
such  as  [19,  24,  25]. 

The  Web  page  will  give  instructions  on  how  to  download  the  software.  Both  compiled 
stand-alone  versions  and  MatLab  scrips  (as  SpecLab  is  now  configured)  will  be  made 
available.  Both  versions  can  be  used  in  Windows  or  Linux-based  environments.  User 
groups  can  be  created  for  the  exchange  of  ideas  and  techniques  involving  SpecLab  and 
similar  software  tools,  data  sets  etc. 

•  SpecLab  manual.  From  the  Web  page  the  user  will  be  able  to  download  a  manual 
that  contains  the  information  on  the  Web  page  itself,  that  is,  a  description  of  all 
subroutines  in  SpecLab,  information  about  downloading  SpecLab,  starting  SpecLab 
etc.  The  manual  will  also  contain  detailed  documentation  of  the  examples  on  the  Web 
page,  as  well  as  detailed  information  about  the  functions  in  SpecLab  and  the  software 
architecture. 

•  Workouts.  The  downloaded  software  will  contain  a  set  of  interactive  examples  that 
will  provide  “guided  tours”  through  specific  applications  of  SpecLab.  The  Phase  I 
version  of  SpecLab  does  this  for  an  aerothermal  data  set.  These  workouts  will  provide 
a  complete  analysis  of  important  data  sets  and  serve  as  an  exended  introduction  to 
SpecLab.  We  envision  at  least  four  workouts:  one  based  on  the  aerothermal  data  set, 
one  based  on  Internet  traffic  data,  and  two  based  on  numerically  generated  data  sets. 

•  Online  help.  At  every  stage  in  the  execution  of  SpecLab  the  user  will  have  access  to 
online  help  in  the  form  of  hyperlinks  to  the  relevant  parts  of  the  formal  documentation. 
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3.2.2  Portability 

Users  will  be  able  to  download  the  package  will  be  available  via  anonymous  ftp  or  directly 
from  the  Web  page.  There  will  be  both  Linux  and  Windows  versions  of  SpecLab.  The  user 
can  choose  to  download  either  a  compiled  version  or  a  version  based  on  Matlab  scripts. 
The  compiled  version  will  not  require  the  user  to  have  Matlab.  The  Matlab  version ,  on  the 
other  hand,  will  enable  the  user  to  tailor  the  Matlab  scripts  to  specialized  processing  tasks. 
This  activity  will  include  replacing  Matlab’s  toolbox  procedures  for  computing  the  wavelet 
transform  with  the  corresponding  routines  in  the  freely  available  software  package  Wavelab 
[13]. 

3.2.3  GUI,  visualization  and  data  manipulation 

•  The  graphical  user  interface  (GUI)  will  be  improved  so  that  it  better  guides  the  user 
through  different  processing  options. 

•  Interactive  parameter  input  will  be  improved.  For  instance,  the  user  will  be  able  to 
select  the  magnitude  of  the  parameters  on  a  continuous  scale  by  using  dials  and  bars 
rather  than  from  the  set  of  discrete  values  in  the  current  version. 

•  The  user  will  be  able  to  employ  the  mouse  to  select  a  subset  of  the  signal  from  a  time 
frequency  plot  and  continue  processing  with  this  subset. 

•  The  routines  for  data  input  and  validation  will  be  improved  as  a  way  to  facilitate  input 
of  data  in  various  formats. 

3.2.4  Adaptation 

Estimation  of  the  locally  stationary  spectral  parameters  is  based  on  a  prior  segmentation  of 
the  data.  This  is  a  key  step.  A  tool  that  supports  the  choice  of  the  prior  segmentation  will  be 
of  great  interest,  not  only  in  SpecLab,  but  also  in  many  other  signal  processing  applications. 
In  version  1  of  SpecLab  the  segmentation  is  done  manually.  In  the  final  version  procedures 
will  be  implemented  to  specify  a  range  of  segmentations.  A  new  component  of  SpecLab  will 
then  look  for  an  optimal  segmentation  within  this  range.  This  new  component  of  SpecLab 
will  in  part  be  based  on  the  framework  set  forth  by  Mallat  et  al.  in  [17].  The  segmentation 
problem  is  a  fundamental  issue  in  signal  processing.  However,  how  to  do  it  well  is  still 
largely  an  open  question.  The  promising  approach  introduced  by  Mallat  et  al.  in  [17]  is 
based  on  an  optimal  compression  of  the  covariance  operator  in  a  suitable  basis.  SpecLab 
provides  a  framework  for  solving  this  important  problem. 

3.2.5  Simulation 

Some  users  may  be  interested  mainly  in  using  SpecLab  for  simulation  of  local  fractional 
Brownian  motion.  Thus,  simulation  of  stochastic  processes  is  another  module  within  SpecLab 
that  will  be  important.  The  current  simulation  algorithm  is  based  on  a  fast  dyadic  refine¬ 
ment  algorithm  [23].  The  proposed  activity  involves  the  following  tasks: 

•  The  simulation  algorithm  will  be  optimized  for  speed.  This  involves  fine  tuning  of  the 
default  parameters  chosen  in  the  simulation  algorithm,  and  also  further  algorithmic 
development.  The  refinement  scheme  in  the  current  simulation  algorithm  will  be 
embedded  in  a  divide-and-conquer  framework  that  will  further  increase  the  efficiency 
of  the  algorithm.  Efficiency  will  be  improved  both  with  respect  to  computational  cost 
and  in  terms  of  use  of  the  memory  hierarchy. 
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•  Rather  than  having  to  use  default  parameters,  the  user  will  be  given  the  option  of 
specifying  a  desired  precision.  With  lower  precision,  i.e.,  larger  deviations  from  the 
specified  power  law,  the  simulation  will  be  faster. 

•  A  separate  menu-driven  tool  that  facilitates  the  simulation  will  be  developed.  The 
user  will  be  able  to  simulate  multiple  realizations  in  parallel  and  carry  out  simulations 
in  a  batch  mode.  This  module  will  make  use  of  the  improved  facilities  for  parameter 
input, 

•  The  user  will  be  given  some  flexibility  in  terms  of  choosing  the  stochastic  model.  It 
will  be  possible  to  simulate  point  samples  of  fractional  Brownian  motion  or  fractional 
Brownian  motion  that  has  been  filtered  so  that  the  simulated  values  correspond  to 
local  integrals  of  pure  fractional  Brownian  motion.  In  addition  the  user  will  be  able 
to  simulate  local  fractional  Brownian  motion  where  power  law  parameters  are  non¬ 
constant. 

•  The  procedure  for  simulation  of  local  fractional  Brownian  motion  will  be  improved.  In 
Version  I  of  SpecLab  this  was  done  with  a  strong  Markov  approximation.  This  will  be 
relaxed  in  version  II  by  means  of  the  divide-and-conquer  step  in  the  new  algorithm. 

3.2.6  Variogram  estimation 

An  important  task  in  applied  signal  processing  is  the  choice  of  an  appropriate  model  for 
the  variogram  or  structure  function.  In  SpecLab  the  variogram  must  be  estimated  for  the 
parameter  processes  that  come  from  the  initial  segmentation.  This  subtask  in  SpecLab  will 
be  developed  further.  The  user  will  be  able  to  modify  parameters  in  a  continuous  fashion, 
exclude  data  points,  modify  the  set  of  points  used  in  the  estimation,  etc.  The  objective  is 
to  provide  a  state-of-the-art  algorithm  for  this  purpose  that  has  at  least  the  functionality  of 
existing  geostatistical  packages.  Since  this  is  an  important  task,  the  user  will  be  given  the 
opportunity  to  carry  out  variogram  estimation  for  an  arbitrarily  specified  data  set,  not  only 
the  parameter  processes. 

3.2.7  Bias,  precision  and  wavelets 

The  current  version  of  SpecLab  is  based  on  Haar  wavelets.  In  the  final  version  the  user  will 
be  given  the  option  of  using  other  types  of  wavelets.  Using  other  wavelets  is  of  particular 
importance  for  extended  signals  where  the  non-stationary  aspect  is  relatively  low  and  the 
desired  precision  for  the  parameter  estimates  is  very  high. 

For  small  sample  sizes  (say  with  only  a  few  spectral  samples)  or  non  Haar  wavelet  based 
estimation  the  bias  in  the  parameter  estimation  must  be  computed  numerically.  This  will  be 
done  a  priori  and  the  bias  looked  up  in  a  table-lookup  procedure.  In  Version  1  of  SpecLab 
an  analytic  expression  for  the  bias  was  used;  this  assumes,  however,  that  the  estimation  is 
based  on  Haar  wavelets. 

In  the  new  version  of  SpecLab  the  user  will  also  be  given  the  option  to  plot  confidence 
bands  for  the  estimated  power  law  parameters.  The  calculation  of  the  variability  of  the 
estimates  will  be  based  on  a  similar  table-lookup  procedure  as  the  one  pertaining  to  the 
bias  correction. 
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3.2.8  Extension  to  exponential  model 

So  far  SpecLab  has  been  tailored  to  local  fractional  Brownian  motion.  In  the  new  version  the 
user  will  be  given  the  option  of  choosing  as  a  model  an  exponentially  correlated  Gaussian 
process  as  well.  This  will  significantly  increase  the  set  of  applications  in  which  SpecLab 
can  be  used  as  a  tool.  The  three  parameters  in  the  exponential  model  (variance,  correlation 
range,  and  white  noise  component)  will  be  modeled  as  being  locally  stationary  and  estimated 
similarly  to  the  way  the  local  power  law  parameters  are  estimated.  This  work  includes 
further  theoretical  development  as  well  as  new  implementation  and  testing  on  some  real 
data  sets.  These  data  sets  can,  for  instance  be  taken  from  the  atmospheric  data  collection 
effort  [34]. 

3.2.9  Incorporating  prior  information 

The  most  important  turbulence  parameter,  the  Hurst  exponent,  takes  on  values  in  the 
range  in  between  zero  and  one.  This  constraint  will  be  included  via  a  non-informative 
prior  distribution  in  a  Bayesian  estimation  scheme.  It  will  also  be  considered  to  include 
more  general  prior  distributions  on  both  turbulence  parameters,  the  Hurst  exponent  and  the 
intensity,  for  better  handling  of  small  data  sets  and  also  for  fast  adaptation. 
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Chapter  4 

Theory 


4.1  Introduction 

Stochastic  processes  that  are  approximately  stationary  and  have  approximately  power  law 
spectral  densities  arise  frequently  in  modeling  atmospheric  turbulence,  financial  data,  geo¬ 
physical  data,  etc.  How  can  we  estimate  the  variable  power  law  behavior  of  the  spectral 
densities  from  data?  The  analysis  will  depend  on  how  we  segment  the  data  and  on  how  we 
choose  the  range  of  scales,  or  frequencies,  over  which  we  look  for  a  power  law  fit.  In  this 
paper  we  address  these  issues  using  fractional  Brownian  motion  as  the  underlying  stochastic 
model  whose  parameters  are  estimated  locally  by  wavelet  scale  spectra.  We  then  apply  the 
theory  to  atmospheric  turbulence  data.  This  data  was  first  analyzed  in  Washburn  et  al.  [34] 
and  subsequently  in  Papanicolaou,  Spina,  and  Washburn  [25];  Rossi,  Kaiser,  and  Washburn 
[30];  Papanicolaou,  Spina,  Rino,  and  Kruger  [19].  We  provide  here  a  general  framework  for 
estimating  local  power  law  processes  with  wavelets  and  a  large  part  of  the  mathematical 
background  needed  for  its  justification. 

In  Section  2  we  introduce  Fractional  Brownian  motion  as  a  model  for  turbulence  and 
discuss  the  wavelet  based  scale  spectrum  that  can  be  used  for  spectral  estimation  of  such 
processes. 

In  order  to  deal  with  approximate  stationarity  and  approximate  power  law  spectra  we 
need  a  good  understanding  of  the  estimation  issues  for  stationary,  power  law  processes, 
like  fractional  Brownian  motion.  Previous  work  on  the  statistics  of  wavelet  scale  spectra 
of  fractional  Brownian  motion  can  be  found  in  Flandrin  [7];  Tewik  and  Kim  [33];  Abry, 
Goncalves  and  Flandrin  [1];  Wornell  [36];  Peltier  and  Levy-Vehel  [27];  Wang,  Cavanaugh 
and  Song  [35].  After  analyzing  the  data,  we  provide  in  Section  6  a  brief  but  complete 
analysis  of  these  statistics  that  is  based  on  a  general  formula  for  the  covariance  of  the 
Haar  wavelet  coefficients  (equation  (4.13)).  A  central  limit  theorem  for  the  estimators  of 
the  power  law  parameters  also  follows  readily  from  the  covariance  formula.  When  Fourier 
spectra  are  used  the  statistical  analysis  of  global  power  law  processes  is  given  in  Robinson 
[32]  and  an  analysis  of  estimators  for  the  exponent  of  the  power  law  is  given  in  Hall,  Koul, 
and  Turlach  [11];  Kent  and  Wood  [16].  We  use  wavelet  scale  spectra  because  they  provide  a 
time-scale  decomposition  of  the  data  that  is  well  suited  to  power  law  processes,  whether  they 
are  stationary  (have  stationary  increments)  or  not.  A  comparison  of  power  law  estimation 
based  on  wavelet  scale  spectra  and  on  Fourier  spectra,  in  the  stationary  case,  is  given  in 
Abry  el  al.  [1]. 

In  Section  3  we  introduce  the  atmospheric  data,  and  carry  out  a  preliminary  scale  spectral 
analysis  using  the  Haar-wavelet  basis.  This  analysis  suggests  that  nonstationary  effects  are 
indeed  important  for  this  data  set.  The  range  of  scales  over  which  the  spectrum  can  be 
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modeled  as  a  power  law,  the  inertial  range,  as  it  is  called  in  turbulence  theory,  varies  with 
spatial  location  and,  in  addition,  the  estimated  power  law  parameters  depend  on  location. 

In  Section  4.4.2  we  address  the  issue  of  choosing  the  range  over  which  to  fit  the  power 
law  by  regression.  We  do  this  by  using  the  fractional  Brownian  motion  as  a  model  locally. 
Since  the  process  is  only  locally  stationary,  power  law  parameters  must  be  estimated  based 
on  relatively  short  spatial  segments.  In  Section  4.4.3  we  show  how  to  remove  the  variability 
of  the  estimated  power  law  parameters  that  is  due  to  the  finiteness  of  the  segments.  The 
filtering  that  we  do  here  is  used  frequently  in  geostatistics  (Ripley  [32];  Cressie  [3]). 

In  Section  5  we  return  to  the  atmospheric  data  set  and  use  the  framework  introduced 
in  Section  4  for  estimation  of  its  scale  spectrum.  We  segment  the  data  into  intervals  over 
which  they  are  approximately  stationary  without  resorting  to  a  full  segmentation  search  as 
in  the  method  of  Mallat,  Papanicolaou  and  Zhang  [17]  that  is  based  on  the  local  cosine 
transform.  To  avoid  the  global  search  we  must  have  a  rough  estimate  of  the  size  of  the 
intervals  of  stationarity,  which  for  the  aerothermal  data  set  we  get  from  a  variogram  analysis 
of  the  wavelet  coefficients  [20].  Another  case  where  the  intervals  of  stationarity  are  known 
approximately  is  analyzed  in  Asch  el  al.  [2]  using  the  local  Fourier  transform. 

Finally  we  show  some  simulations  from  the  estimated  model  that  assess  the  overall 
relevance  of  our  analysis.  The  main  point  of  our  analysis  is  that  we  are  able  to  identify  the 
local  variations  of  the  power  law  parameters,  the  Kolmogorov  turbulence  law,  which  arise 
from  large  scale  atmospheric  phenomena.  The  analysis  should  be  applicable  to  other  data 
sets,  financial  data  for  example,  where  departure  from  stationarity  needs  to  be  quantified. 


4.2  Scale  spectrum  of  Fractional  Brownian  motion 

4.2.1  Fractional  Brownian  motion 

We  shall  model  ‘pure’  power  law  processes  by  fractional  Brownian  motion  (fBm),  {Bh{x)\  x  > 
0},  introduced  by  Mandelbrot  and  Van-Ness  [20].  It  is  a  Gaussian  process  with  mean  zero, 
stationary  increments  and  covariance 

E[BH{x)BH{y)}  =  y(M2"  +  \y\2H-\x  -y\2H)  (4.1) 

with  0  <  H  <  1  and  a  parameters.  Its  structure  function  is 

E[{BH(x)  -  BH(x  -  Ax))2]  =  a2\Ax\2H, 

and  it  is  conditioned  to  be  zero  at  the  origin:  Bh(0)  =  0.  The  so  called  Hurst  exponent 
H  determines  the  correlation  of  the  increments.  The  covariance  of  future  increments  with 
past  ones  is 


Ph(&x)  =  E[{Bh{x)  -  Bh(x  -  Ax)) 

x{Bh(x  +  Ax)-Bh(x))]  =  <r2(  22"-1-1)|Ax|2", 

which  is  independent  of  x.  When  H  >  1/2  this  quantity  is  positive  so  if  the  past  incre¬ 
ment  is  positive,  then  on  average  the  future  increment  will  be  positive.  Feder  [5]  calls  this 
persistence.  When  H  <  1/2  we  have  an  antipersistent  process  with  a  positive  increment  in 
the  past  making  a  positive  increment  in  the  future  less  likely.  Ordinary  Brownian  motion 
corresponds  to  H  =  1/2.  In  this  case  future  and  past  increments  are  independent.  Of 
special  interest  here  is  H  —  1/3  which  corresponds  to  the  Kolmogorov  scaling  law  (Monin 
and  Yaglom  [21]). 
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Fractional  Brownian  motion  has  stationary  increments  but  is  not  itself  stationary.  How¬ 
ever,  as  shown  in  Abry  et  al.  [1]  for  example,  it  is  possible  to  assign  a  pseudo-spectrum  to 
it  by  cutting  off  the  low  frequencies.  Let 

X  =  Bh  *  V* 

with  star  denoting  convolution  and  ib  a  function  that  integrates  to  zero  so  that  its  Fourier 
transform  vanishes  at  zero  frequency.  The  process  X  is  stationary  and  its  power  spectrum 
is 

Px  <X  <72|/|-(2if+1)|^(/)|2. 

Since  we  usually  observe  power  law  processes  through  a  filter  that  cuts  off  very  low  frequen¬ 
cies,  we  can  associate  with  Bh  the  power  law  spectrum 

PBh  oc  a2|/|-<2//+1>  . 

In  the  Kolmogorov  case  H  —  1/3  the  spectrum  is  <r2|/|”5/3  over  some  range  of  frequencies, 
called  the  inertial  range. 

Fractional  Brownian  motion  is  self-similar  since  Bh{%)  and  aH  Bh(x/ a)  have  the  same 
finite  dimensional  distributions  for  all  a.  We  next  discuss  estimation  of  fBm,  or  a  pure  power 
law  process.  In  Section  4.4  we  generalize  the  estimation  procedure  to  locally  stationary 
power  law  processes. 


4.2.2  Haar  wavelets  and  scale  spectrum 


We  want  to  carry  out  a  spectral  analysis  of  the  process,  given  as  a  finite  set  of  data,  and  to  fit 
the  estimated  spectrum  to  a  power  law.  For  this  purpose  scale  spectra,  rather  than  Fourier 
spectra,  will  be  used.  Scale  spectra  are  a  natural  and  flexible  tool  for  self-similar  processes 
(Kaiser  [15];  Abry  et  al.  [1]).  The  scale  spectrum  is  defined  in  terms  of  the  coefficients  of  the 
data  in  a  wavelet  basis.  We  will  use  the  Haar  wavelet  basis  although  other  bases  could  have 
been  used  as  well  (Goncalves  and  Abry  [9]).  Haar  wavelet  based  estimators  are  sometimes 
related  to  classical  ones,  as  is  in  particular  the  Allan  variance  estimator  for  the  slope  of  the 
log  of  the  spectral  density  (Percival  and  Guttorp  [28]). 

Let  Y  denote  the  process  for  which  we  want  to  compute  the  Haar  wavelet  coefficients. 
Denote  the  approximation  coefficients  at  level  zero  by  X  =  (a0(l),  ao(2), ...,  ao(2M)).  These 
are  defined  by 

flX 

a0(n)  =  /  Y{x)dx  . 

Jn-l 

Given  this  level  zero  representation,  the  data,  we  construct  successively  its  wavelet  coeffi¬ 
cients  with  respect  to  the  Haar  basis  as  follows.  Let 


ax(n) 


d\  (n) 


V2 


(a0(2n)  +  a0(2n  -  1)) 


—  (a0(2n)  -  a0(2n  -  1)) 


for  n=l,  2...,  2M~l 


(4.2) 


be  the  smoothed  signal  and  its  fluctuation,  or  detail,  at  the  finest  scale.  Note  that  the  detail 
vector  d\  contains  every  other  successive  difference  of  the  data.  This  process  of  averaging 
and  differencing  can  be  continued  by  defining 


ao(n)  = 


d‘2(n)  = 


~{ai{2n)  +  ai(2 n  -  1)) 
(ai(2n)  -  ai(2n  -  1)) 


for  ri  =  1.2...,2m"2 
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and  in  general 


aj  (n) 
dj(n) 


-j=(aj-i{2n)  +  Oj_i(2n-  1)) 
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[cij-i  (2n)  —  aj_i  (2n  —  1))  ,  for  n  =  l,2...,2 


M-j 


for  j  =  1  The  data  vector  X  can  then  be  reconstructed  from 

since  from  equations  (4.2)  we  have 

ao(2n)  =  -^=(ai(n)  +di(n)) 

ao(2n-l)  =  -i=(ai(n)  -  c?i(n))  ,  for  n  —  1,  2...,  2M_1 

v  2 

and  now  ai  can  be  replaced  by  sums  and  differences  of  a*?  and  d2>  etc. 

The  detail  coefficients  at  level  j  can  alternatively  be  expressed  as 

i  r°° 

dj(n)  =  —j=  /  ip(x/23  -  n)Y(x)dx 

v2J  7—00 

with  the  mother  wavelet  defined  by 

-1  if  -1  <x  <  -1/2 
ib(x)  =  1  if  —  1/2  <  x  <  0 

0  otherwise 


The  difference  coefficients  correspond  to  probing  the  process  at  different  scales  and  locations, 
with  n  representing  location  and  j  scale.  From  the  self-similarity  of  fractional  Brownian 
motion  it  follows  that  for  such  a  process  E[dj(n )2]  oc 

The  scale  spectrum  of  A"  relative  to  the  Haar  wavelet  basis  is  the  sequence  Sj  defined 
by 


Sj  =  2A7 -  E  (d»)2  -  3  =  1, 2, M.  (4.3) 

rt—  1 

For  fractional  Brownian  motion  the  log  of  the  scale  spectrum  is  approximately  linear  in  the 
scale  j.  As  discussed  further  in  Section  6  this  can  be  used  to  estimate  the  parameters  of  the 
process  by  regression  It  is  easily  verified  from  the  definitions  above  that  the  l 2  norm  of  the 
data  vector  X  can  be  written  as 

2m  m 

X>0 (n))2  =  (aM)2  +  E2M_JA  . 

n— 1  j= 1 

which  is  a  way  of  expressing  the  orthogonality  of  the  decomposition  of  X  into  gm  and  the 
dj,  j  =  1, M. 

The  scale  spectral  point  Sj  is  the  mean  square  of  the  detail  coefficients  at  scale  j.  The 
spectrum  can  therefore  be  interpreted  as  the  energy  of  the  signal  in  the  different  scales 
relative  to  the  Haar  wavelet  basis.  Consider  data  containing  information  only  at  the  finest 
scale:  X  =  {1,  —  1,  1,  —1,  •— }.  Then  only  the  d\{n)  coefficients  are  non-zero.  Hence, 

Sj  roughly  corresponds  to  the  energy  at  times  the  finest  scale. 

Note  that  the  spatial  support  of  the  integrals  of  the  process  defining  the  wavelet  coef¬ 
ficients  at  a  certain  scale  is  adapted  to  the  particular  scale.  The  short  scales,  that  provide 
high  frequency  information,  are  defined  in  terms  of  consecutive  integrals  of  narrow  support. 
A  plot  of  these  difference  coefficients  reveals  information  about  how  the  high  frequency 
content  of  the  data  change  with  location. 
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4.3  The  aerothermal  data 

4.3.1  Temperature  and  index  of  refraction  fluctuations 

Our  objective  is  to  estimate  the  spectrum  of  temperature  data  taken  from  the  upper  atmo¬ 
sphere.  One  motivation  for  our  modeling  and  estimation  is  to  be  able  to  generate  synthetic 
random  media  for  the  simulation  of  wave  propagation  in  a  turbulent  atmosphere.  Before 
we  analyze  the  aerothermal  data  we  explain  briefly  how  velocity,  temperature  and  index 
of  refraction  fluctuations  are  related  to  each  other  and  how  the  Kolmogorov  scaling  theory 
enters  in  their  description  (Monin  and  Yaglom  [21];  Strohbehn  [31]). 

The  Kolomogorov  theory  is  a  phenomenological  statistical  description  of  the  velocity  field 
in  the  atmosphere.  Based  on  a  scaling  argument,  the  mean-square  velocity  differences  are 
described  in  a  universal  manner  over  a  rather  broad  range  of  spatial  scales,  the  inertial  range. 

If  we  assume  homogeneity,  isotropy  and  incompressibility  the  result  is  that  the  structure 
function  of  the  velocity  has  the  form 

£[(vx(x0+x)-vx(x0))2]  =  Cl  |x|2/3  (4.4) 

Here  vx  is  the  velocity  in  the  direction  of  the  displacement  x.  The  parameter  C'l  is  the 
velocity  structure  constant,  a  measure  of  the  amount  of  energy  in  the  inertial  range,  which 
is  typically  confined  between  an  inner  scale  Iq  and  an  outer  scale  To. 

The  connection  to  temperature  is  through  the  theory  of  convection-diffusion  of  a  passive 
scalar.  It  turns  out  (Monin- Yaglom  [21])  that  temperature  fluctuations  have  also  a  structure 
function  with  Kolmogorov  2/3  scaling. 

At  optical  frequencies  the  variations  in  the  index  of  refraction  Sn  are  approximately 

Cqn 

Sn  =  -79 P—r  x  1CT6 

with  P  atmospheric  pressure  in  millibars  and  T  temperature  in  degrees  Kelvin.  Thus,  in 
order  to  describe  the  spectrum  of  fluctuations  in  the  refractive  index  we  need  only  the 
spectrum  of  the  fluctuations  in  the  temperature.  Wavelet  scale  spectra  are  used  to  analyze 
turbulence  data  in  Hudgins,  Friehe  and  Mayer  [14]. 

The  above  model  for  refractive  index  fluctuations  is  used  extensively  to  model  atmo¬ 
spheric  wave  propagation.  Variations  in  the  turbulent  character  of  the  medium  are  typically 
modeled  through  variations  in  only  and  not  in  the  exponent,  which  is  fixed  at  2/3. 

4.3.2  Data  analysis  of  aerothermal  data 

We  will  analyze  temperature  data  obtained  by  the  Air  Force  high-altitude  laser  propagation 
and  turbulence  data  collection  effort.  For  a  detailed  discussion  of  recording  procedures 
and  analysis  see  Washburn  et  al.  [34] ;  Papanicolaou  et  al.  [34];  Rossi  et  al.  [30].  A 
strong  effort  was  made  to  provide  high  quality  data  that  could  be  used  to  characterize  the 
turbulent  atmosphere.  This  unique  data  set  cannot  be  maid  generally  available,  but  the 
above  references  provide  together  a  fairly  detailed  account  for  the  data.  Here  we  take  the 
data  as  our  starting  point  and  examine  what  analysis  reveals  about  their  structure.  We  are 
interested  in  particular  in  accounting  properly  for  nonstationary  effects.  All  calculations  are 
carried  out  in  MATLAB  on  a  Silicon  Graphics  workstation. 

In  Figure  4.1  we  show  the  temperature  data,  which  has  approximately  4.2  x  106  points. 

The  spatial  resolution  is  approximately  2cm.  The  data  is,  of  course,  quite  noisy  and  it  is 
part  of  our  task  to  remove  spurious  noise  effects  in  the  estimation  process. 

In  Figure  4.2  we  show  the  first  21  km  of  the  temperature  data  along  with  di,  d5,  d7,  d9,  dn,  an- 

The  detail  coefficients  dj  carry  information  about  the  data  on  larger  scales  as  j  increases.  For 
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Data;  subsampled  1/50 


Figure  4.1:  Top  figure:  the  raw  temperature  data.  The  spatial  resolution  of  the  data  is 
approximately  2cm. 


example,  dj  shows  successive  differences  of  the  temperature  over  distances  of  2.56m,  after 
averaging  over  successive  segments  of  length  1.28.  Thus,  as  j  increases  the  data  are  lowpass 
filtered  with  a  filter  of  decreasing  bandwidth  and  then  subsampled  before  the  differences  are 
formed.  It  is  clear  that  the  visible  periodic  components  seen  in  the  d\  coefficients  cannot  be 
attributed  to  turbulence  or  larger  coherent  structures  in  the  atmosphere.  In  the  estimation 
of  atmospheric  turbulence  parameters  these  spurious  features  must  be  suppressed.  In  c?5,  •  •  • 
the  high  frequency  periodic  noise  component  seen  above  has  been  effectively  suppressed  by 
the  lowpass  filtering. 

Our  main  objective  is  now  to  examine  whether  the  data  can  be  modeled  well  by  a  ‘power 
law’  model  over  a  subrange  of  scales.  Recall  that  for  such  a  process,  the  scale  spectrum  Sj 
is  linear  in  a  log-log  plot,  when  the  record  is  very  long.  In  Figure  4.3  we  show  log-log 
plots  of  scale  spectra  over  nonoverlapping  segments  of  the  data,  with  each  segment  having 
length  1.3km.  The  top  plot  corresponds  to  the  first  half  of  the  data  whereas  the  bottom 
plot  corresponds  to  the  second  half.  Each  plot  contains  scale  spectral  points  over  16  scales, 
1  <  j <  16.  This  corresponds  to  length  scales  from  li  =  .04m  to  /i6  =  1.3A;m.  We  take  as 
abscissa  the  spatial  frequency  Kj 

Kj  =  [rad/m]  =  1007t2“j  [rad/m].  (4.5) 

Note  that  the  scale  spectra  show  a  distinct  departure  from  power  law  behavior  for  scales 
below  one  meter  ( K  =  6).  For  larger  scales,  power  law  behavior  may  be  considered  in  the 
range  between  2.5m  to  80m  (.08  <  K  <  2.5)  approximately,  which  corresponds  to  the  detail 
coefficients  dj  to  di2-  The  departure  from  power  law  behavior  for  the  shortest  scales  is 
partly  due  to  the  influence  of  measurement  noise.  The  first  half  of  the  data  is  less  energetic 
than  the  second  half  and  hence  more  noisy.  The  estimated  intercept  and  slope  for  the  log 
scale  spectra  depend  on  the  particular  segment  of  data  used.  We  want  to  identify  the  part 
of  this  variability  that  is  due  to  the  nonstationarity  of  the  process  and  minimize  variability 
due  to  estimation  errors. 

For  fractional  Brownian  motion  the  wavelet  coefficients  will  be  normally  distributed. 
In  Figure  4.4  we  show  a  histogram  of  the  wavelet  coefficients  at  scale  8,  normalized  by  a 
local  estimate  of  the  variance.  The  dashed  line  corresponds  to  a  Gaussian  distribution.  In 
the  bottom  plot  we  show  the  autocorrelation  of  these  wavelet  coefficients,  normalized  by 
their  variance.  We  plot  the  correlation  in  terms  of  the  empirical  variogram  computed,  as 
in  (4.11)  below.  The  dashed  line  is  the  theoretical  correlation  for  fBm  with  Kolmogorov 
scaling,  H  =  1/3  (as  defined  in  (4.13)).  In  Figure  4.5  we  show  an  estimate  of  the  structure 
function  for  log (dj(n)2);  7  <  j  <  9.  If  the  process  had  a  pure  power  law  spectrum  over  these 
scales  we  would  see  only  a  very  small  correlation.  However,  these  quantities  are  correlated 
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Lovel  of  detail  coefficients 


Figure  4.2:  Haar  wavelet  coefficients  for  the  first  21km  of  the  data.  The  top  left  figure  is 
the  temperature  data.  The  one  below  it  is  the  d\  detail  coefficients.  The  third  from  the  top 
is  the  g?3  and  the  bottom  the  d5  detail  coefficients.  Note  that  the  4km  burst  seen  in  di  is 
not  visible  any  more  at  the  d$  level,  after  four  successive  averagings  of  the  data.  The  top 
figure  on  the  right  is  the  an  coefficients  and  below'  it  dn,d§,d7. 


on  the  order  of  km ,  indicating  that  we  have  a  local  power  law  structure. 

Before  we  continue  with  the  scale  spectral  analysis  we  note  the  following. 

•  Noise  bursts  in  the  data  enter  only  in  some  of  the  detail  coefficients,  as  can  be  seen 
from  Figure  4.2.  The  criterion  for  selecting  of  the  inertial  range,  implemented  in  the 
next  section,  will  automatically  restrict  the  scale  range. 

•  The  temperature  data  are  not  a  stationary  time  series  and  they  do  not  have  stationary 
increments.  Computing  scale  spectra  over  long  segments,  in  which  the  process  cannot 
be  taken  as  stationary,  gives  a  quantity  that  is  hard  to  interpret.  Even  though  the 
average  slope  of  the  log  scale  spectra  over  several  segments  is  close  to  —5/3,  as  the 
theory  of  turbulence  predicts,  there  is  a  lot  of  variability.  A  local  power  law  model 
of  the  kind  discussed  in  Section  4.4  is  likely  to  fit  the  data  better  than  the  idealized 
power  law  model,  with  stationary  increments. 

We  will  estimate  the  scale  spectrum  of  the  data  when  we  model  it  as  a  local  power 
law,  using  the  method  described  in  Section  4  To  start  the  estimation  we  must  have  rough 
estimates  for  the  intervals  of  stationarity  and  for  the  exponent  of  the  pow'er  law.  In  view  of 
Figures  4.3,  4.4  and  4.5  we  choose  these  estimates  to  be 

•  (f)-1  =  2  km~l 

•  Hq  =  1/3  (Kolmogorov  scaling)  . 

In  the  next  section  we  develop  a  framework  for  estimation  of  local  power  law  processes  with 
varying  inertial  range  and  power  law  parameters.  In  Sections  4.5.1  and  4.5.2  we  carry  out 
the  estimation.  First  we  estimate  the  inertial  range  and  then  the  power  law  parameters. 
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Scalo  Spectra 


Figure  4.3:  Scale  spectra  of  1.3A;m  nonoverlapping  segments  of  the  data  (216  points  per 
segment)  obtained  from  the  Haar  wavelet  decomposition.  The  top  plot  corresponds  to  the 
first  half  of  the  data  whereas  the  bottom  plot  corresponds  to  the  second  half.  The  dashed 
line  has  slope  —5/3  as  in  Kolmogorov  spectra.  The  solid  line  is  the  average  over  the  scale 
spectra  of  the  different  segments.  After  the  averaging  is  done  the  scale  spectra  are  plotted 
in  log- log  format. 

4.4  Estimation  of  local  power  law  processes 

4.4.1  Modeling  and  segmentation 

The  fractional  Brownian  motion  described  in  Section  2  is  an  idealization  of  a  process  with 
power  law  spectrum.  It  has  stationary  increments  and  if  the  power  law  is  truncated  at 
low  frequencies  then  it  is  itself  stationary.  In  most  physical  or  financial  applications  where 
power  law  behavior  is  expected  locally,  the  power  law  parameters  will  vary  so  the  process 
cannot  be  stationary  in  the  large.  The  aerothermal  data  provide  an  example.  We  would 
like  to  estimate  the  power  law  parameters  over  segments  that  are  short  enough  so  that  they 
can  be  taken  as  constant  but  long  enough  so  that  their  statistical  estimates  are  stable.  How 
are  we  to  decide  what  is  a  good  segmentation  of  the  data  in  this  vague  sense?  This  is  a 
very  difficult  problem  that  is  rarely  addressed  in  theoretical  or  applied  studies  of  spectral 
estimation.  In  Mallat  et  al.  [17]  a  class  of  locally  stationary  processes  is  introduced  and 
an  algorithm  for  finding  intervals  of  approximate  stationarity  is  developed.  This  method 
is  based  on  an  exhaustive  search  for  an  optimal  segmentation,  in  a  suitable  sense,  into 
intervals  of  stationarity.  It  does  not  work  so  well  when  the  intervals  of  stationarity  are  all 
of  roughly  the  same  size,  and  when  dealing  with  turbulence  data.  Thus,  here  we  want  to 
follow  a  somewhat  different  approach  where  we  can  take  advantage  of  prior  information 
about  intervals  of  stationarity,  as  is  the  case  of  the  aerothermal  data. 

To  fix  ideas  we  will  consider  estimation  of  parameters  for  a  multifractional  Brownian 
motion  (mBm).  This  is  a  generalization  of  fractional  Browmian  motion  where  the  parameters 
vary  with  location.  In  the  time  domain  they  have  the  representation,  (  Goncalves  and 
Flandrin  [9];  Peltier  and  Levy-Vehel  [27]) 

B,(x)  =  f TiiTXTTlo\  f  [(»-^,(J)-1/2-(-*)",(x)~1/3]^)  (4-6) 

I  (h;(x)  +  1/2)  J_ 30 
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Marginal  of  d8  normalized:  solid  line  is  Normal 


Figure  4.4:  The  top  plot  shows  the  empirical  distribution  of  d$.  The  solid  line  is  a  Gaussian 
distribution.  The  bottom  plot  shows  the  correlation  of  these  wavelet  coefficients.  The 
dashed  line  is  the  theoretical  correlation  for  fBm  with  H  —  1/3. 


+  [X {x  -  s)H^x)~l/2dW(s) 

Jo 

with  W  the  standard  Brownian  motion  and 

Hs(x)  =  H{ex)  (4.7) 

as(x)  =  <r(ex)  (4.8) 

where  H(-)  &  cx(-)  can  be  deterministic  or  random.  For  example,  they  can  be  stationary 
stochastic  processes  with  smooth  paths  and  decaying  correlation  functions,  independent 
of  the  Brownian  path  W(-).  When  H  and  a  are  constants  then  (4.6)  is  a  time  domain 
representation  of  the  usual  fractional  Brownian  motion. 

The  parameter  e~l  is  a  measure  of  the  interval  of  stationarity  in  the  sense  that,  as  in 
Goncalves  and  Flandrin  [9], 

E[(Bs(x)  -  B*(x  —  Ax))2]  «  a?(x)|Ax|2ift'^ 

for  eAx  small.  This  means  that  for  scales  that  are  small  compared  to  the  interval  of 
stationarity,  the  processes  Be{x)  behaves  locally  like  a  fractional  Brownian  motion  with 
parameters  frozen  at  x. 

In  Section  4.4.3  we  will  describe  a  method  that  removes  dependence  of  the  estimated 
parameters  on  a  prior  estimate  of  the  interval  of  stationarity.  The  basic  idea  is  that  we  know 
roughly  what  1  je  should  be  and  how  the  estimation  errors  behave  if  the  underlying  process 
is  fractional  Brownian  motion,  given  the  segment  size.  We  then  filter  out  these  estimation 
errors  and  get  estimates  that  do  not  depend  sensitively  on  the  prior  choice  of  the  size  of  the 
interval  of  stationarity. 

Another  modification  that  is  necessary  in  the  nonstationary  case  is  the  identification  of 
the  inertial  range  over  each  interval  of  approximate  stationarity.  The  error  in  the  power 
law  fit  to  the  scale  spectra,  over  a  segment  of  data,  depends  on  the  scale  range  that  is  used. 
How  do  we  select  a  range  of  scales  for  which  the  error  in  the  fit  is  acceptable?  In  the  next 
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Structure  function  for  log(cfj(n) 2); 


Figure  4.5:  The  figure  shows  an  estimate  of  the  structure  function  of  lo g{dj(n)2)  as  function 
of  nAx  and  with  7  <  j  <  9.  The  magnitude  of  these  quantities  determine  the  scale  spectrum 
and  they  decorrelate  over  a  scale  on  the  order  of  km.  The  dotted,  dashed  and  solid  lines 
corresponds  respectively  to  scales  7,8  and  9. 


section  we  introduce  a  criterion  for  selecting  the  range  based  on  comparison  with  an  ideal 
fractional  Brownian  motion. 

4.4.2  Estimation  of  the  inertial  range 

The  power  law  behavior  that  we  want  to  identify  in  the  data  is  necessarily  restricted  to  a 
finite  set  of  scales,  the  inertial  range,  in  each  segment  of  stationarity. 

Let  2m  denote  the  length  of  the  data  vector  over  the  segment  under  consideration.  For 
an  ideal  fractional  Brownian  motion  with  Hurst  exponent  H  let  S*,  i  —  1, . . . ,  M  be  its 
wavelet  scale  spectrum.  Then,  for  each  subrange  {z‘i,  •  -  * ,  i2}  the  measure  of  misfit 

r(ti,t3)  =  j>g2(Si)-log2(&]a>  (4-9) 

i=ii 

with  Si  the  power  law  estimated  by  weighted  least  squares,  is  a  random  variable  whose  law 
can  be  computed  analytically  in  principle  or  numerically.  It  depends  weakly  on  H  so  we  fix 
it  to  equal  some  rough  estimate  H  =  H0.  For  the  aerothermal  data  of  Section  4.3  we  take 
H0  =  1/3,  the  Kolmogorov  exponent.  Let  R(ix,i 2)  be  the  90th  percentile  of  the  distribution 
of  r{i\ ,  i2)  with  the  value  of  the  Hurst  parameter  equal  to  Ho.  The  scale  range  is  now  chosen 
as  the  largest  range  for  which  r(ii,i2)  <  R(i  1,^2)-  Here  r(ii,i2)  is  the  value  of 

the  error  (4.9)  obtained  from  the  actual  data. 

In  order  to  obtain  stable  power  law  parameter  estimates  in  the  regression  we  need  a  min¬ 
imum  number  of  scale  range  points.  In  the  application  to  the  aerothermal  data  introduced 
in  Section  3  we  choose  to  model  in  terms  of  a  power  law  only  over  segments  for  which  there 
are  5  or  more  points  in  the  estimated  scale  range. 

4.4.3  Segmentation  independent  power  law  estimation 

Given  a  fixed  segmentation  of  the  data  into  approximate  intervals  of  stationarity  that  is 
based  on  some  prior  information,  we  first  calculate  the  inertial  range  as  described  in  the 
previous  section  and  then  do  the  power  law  fit.  The  power  law  parameters  are  obtained  by 
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weighted  linear  regression  of  the  log  scale  spectrum,  see  Section  4.7.4.  The  estimated  power 
law  parameters  depend  on  the  segmentation.  We  will  now  describe  a  method  with  which 
this  dependence  can  be  removed. 

The  idea  is  to  do  a  filtering  of  the  parameter  estimates  in  order  to  remove  the  variability 
that  is  segmentation  dependent.  From  the  theory  of  the  power  law  estimators  for  fractional 
Brownian  motion  we  know  that  the  slope  estimator  has  the  form 

p*  -  pi  +  wi  (4.10) 

where  pl  is  the  slope  for  the  i’th  segment  and  wl  the  fluctuation  due  to  the  finiteness  of 
the  segment.  We  cannot  take  large  segments,  that  reduce  this  error,  because  we  are  limited 
by  the  nonstationarity.  The  errors  wx  are  essentially  uncorrelated  over  different  segments 
and  we  have  to  construct  a  filter  that  will  predict  pl  from  the  estimates  p%  by  removing  the 
noise  w\  to  the  extent  possible.  We  assume  that  p1  is  itself  a  stationary  stochastic  process, 
independent  of  the  Brownian  motion  that  generates  the  fractional  Brownian  motion.  The  pl 
are  the  intrinsic  variation  of  the  power  in  the  locally  stationary  fractional  Brownian  motion. 
The  correlation  length  of  the  pl  must  be  longer  than  the  segments  used  in  estimating  them 
in  order  to  have  approximate  stationarity  relative  to  the  segmentation. 

It  can  be  shown  using  the  results  from  Section  6  that  wl  is  close  to  being  a  white  process. 
This  is  important  because  it  allows  us  to  estimate  the  autocovariance  of  px  from  the  estimates 
pl,  using  the  variogram,  for  example.  Given  estimates  of  the  autocovariance  of  pl  and  the 
variance  of  wl  we  can  design  a  minimum  mean  square  error  predictor  or  smoothing  filter  for 

p\ 

The  minimum  variance  unbiased  filter  for  prediction  of  the  parameter  processes  is  as 
follows.  We  describe  it  in  the  context  of  the  slope  parameter,  but  the  filtering  of  the  log 
intercept  is  completely  analogous.  Let  P  =  (pi)  be  the  vector  of  estimates,  P  =  (pi)  the 
realization  of  the  slope  process  and  P  =  (p)  the  constant  mean.  Then  the  filter  F  is  a  matrix 
that  transforms  P  into  TP  in  such  a  way  that 

£[||rp-p||2] 

is  minimized  over  all  matrices  T  that  also  preserve  the  mean  P  of  P,  that  is  FP  =  P.  Let 
Cp  be  the  covariance  matrix  of  P.  Then  it  easily  follows  that 

r  =  (Cp  +  cj-^  +  u7®:?] 

where  the  vector  u  —  (u^)  is  given  by 

Pi-P  T(Cp  +  Cw)~1CPii 
Ui  -  P  T(Cp  +  Cw)~'P 

and  Cw  is  the  diagonal  covariance  matrix  of  the  estimation  errors  wl .  Here  Cp<i  is  the  i— th 
column  of  the  matrix  Cp  and  the  superscript  T  stands  for  transpose.  The  slope  and  log 
intercept  processes  are  filtered  separately.  Filtering  of  this  kind  is  discussed  in  Ripley  [32] 
and  Cressie  [3]  for  example. 

Since  the  effect  of  the  sample  noise  wl  in  (4.10)  is  largely  removed  by  the  filtering 
procedure,  the  estimates  of  the  parameter  processes  will  be  essentially  independent  of  the 
prior  choice  of  segmentation. 

4.5  Application  to  aerothermal  data 

4.5.1  Estimation  of  inertial  range 

We  estimate  the  set  of  scales  where  the  process  can  be  modeled  as  a  'power  law’.  Figure 
4.3  shows  that  it  varies  with  location.  That  is,  the  set  of  scales  where  the  scale  spectrum  is 
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approximately  linear  is  location  dependent.  There  are  two  main  reasons  for  this.  First,  the 
scale  range  where  the  physical  process  has  power  law  spectrum  varies  and  depends  on  the 
local  intensity  of  the  turbulence.  Second,  the  set  of  scales  that  are  affected  by  measurement 
noise  varies  depending  also  on  the  intensity  of  the  turbulence. 

We  use  the  scheme  described  in  Section  4.4.2  to  estimate  the  inertial  range.  First, 
we  choose  a  segmentation  that  is  short  relative  to  the  prior  estimate  of  the  interval  of 
stationarity.  We  choose  segments  of  length  215  points.  Then  we  apply  the  algorithm  of 
Section  4.4.2.  For  all  scale  ranges  we  measure  the  difference  between  the  scale  spectrum 
computed  from  the  data  and  the  fitted  power  law  and  choose  the  estimated  scale  range  as 
the  largest  one  for  which  this  measure  is  within  the  90th  percentile  of  the  corresponding 
measure  for  a  realization  of  fractional  Brownian  motion. 

The  estimated  ‘effective  inertial  ranges’  are  shown  in  Figure  4.6  by  the  vertical  solid 
lines.  For  each  segment  there  is  one  vertical  line  showing  the  scales  included  in  the  corre¬ 
sponding  effective  inertial  range.  The  segments  consist  of  2lD  points  and  the  maximum  scale 
considered  is  therefore  14.  Note  that  the  more  energetic  second  half  of  the  data  displays 
a  more  consistent  power  law  behavior,  as  one  would  expect  on  physical  grounds  in  fully 
developed  turbulence. 


Scale  ranges  with  residuals  consistent  with  model 


[km} 


Figure  4.6:  The  vertical  lines  in  the  figure  are  the  estimates  of  the  location  dependent  inertial 
ranges.  That  is,  the  scales  over  which  the  observed  process  is  approximately  a  power  law. 
For  each  segment  of  length  2io  points  there  is  one  vertical  line  showing  the  scales  included 
in  the  inertial  range  estimate. 

The  spectra  computed  with  the  estimated  location  dependent  inertial  ranges  are  shown 
in  Figure  4.7,  the  top  plot.  They  are  plotted  in  a  log-log  format  and  shifted  vertically  to 
coincide  at  a  center  scale.  We  include  only  spectra  that  have  at  least  5  points.  Note  the 
fluctuations  in  the  computed  slopes.  In  Figure  4.7,  bottom  plot,  we  show  the  average  of 
the  spectra.  The  averaging  is  carried  out  before  the  log  transformation.  Note  the  excellent 
match  with  the  Kolmogorov  scaling  law  shown  by  the  solid  line.  There  is  a  slight  deviation 
at  the  5th  scale.  This  corresponds  approximately  to  the  onset  of  the  measurement  noise 
caused  by  the  airplane,  see  Figure  4.3. 

4.5.2  Estimation  of  power  law  parameters 

In  the  previous  Section  we  estimated  the  inertial  range  for  the  scale  spectrum,  relative 
to  a  fixed  segmentation.  We  will  now  estimate  the  power  law  parameters  over  this  set  of 
constrained  scales.  We  use  the  method  described  in  Section  4.4.3  for  this  purpose.  Our 
main  objective  is  to  capture  the  intrinsic  variation  in  the  parameters.  If  we  choose  a  coarse 
segmentation  the  variability  of  the  estimates  will  be  small.  However,  the  estimated  power 
law  based  on  the  scale  spectrum  may  in  this  case  be  an  average  of  power  law  parameters 
that  vary  within  the  section.  This  is  illustrated  in  Figure  4.7.  The  average  power  law 
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Figure  4.7:  The  top  plot  shows  the  computed  scale  spectra  in  each  subsection,  dotted  lines. 
The  estimates  are  based  on  the  scale  ranges  shown  in  Figure  4.6.  The  bottom  plot  are  the 
average  of  the  (untransformed)  spectra.  Note  the  good  match  with  the  Kolmogorov  scaling 
law  shown  by  the  solid  line. 


over  the  whole  data  set  is  very  close  to  the  Kolmogorov  scaling  law.  However,  within  each 
section  (see  top  plot)  such  a  model  must  clearly  be  rejected.  Therefore  we  must  choose  a 
segmentation  that  is  not  much  larger  than  the  interval  of  stationarity.  The  estimates  should 
be  independent  of  segmentation.  That  is,  if  we  shorten  the  segments  this  should  not  lead 
to  an  increase  in  the  variability  of  the  estimated  parameters. 

To  show  that  this  is  achieved  by  filtering  we  choose  three  segmentation  lengths  in  addition 
to  the  one  used  in  the  previous  section.  We  choose  the  following  segmentations 

160m  (213  points),  327m  (214),  655m,  (215),  1.31/cm  (216). 

There  are  512,  256,  128,  64  nonoverlapping  segments,  respectively,  in  each  case.  For  each 
segmentation  we  estimate  inertial  ranges  for  the  scale  spectrum  as  above.  In  Figure  4.6  we 
show  the  estimated  ranges  corresponding  to  segments  of  length  655m. 

We  use  the  linear  least  squares  regression 

log  Sj  «  cl  +pilog2(^) 

for  each  segment  z,  with  j  the  scale  in  suitable  units.  In  Section  4.7.4  we  analyze  the 
regression  (4.11).  The  results  of  the  estimation  are  shown  in  Figure  4.8.  If  the  inertial  range 
has  fewer  than  5  spectral  points  we  do  not  estimate  power  law  parameters  and  leave  a  gap 
in  the  figure. 

We  see  from  the  figure  that  the  estimated  slopes,  p,  and  the  log  intercepts,  c,  vary  consid¬ 
erably  over  the  80  kilometer  data  set.  They  also  depend  on  the  segmentation,  with  the  finer 
ones  having  larger  fluctuations  (dotted  lines).  Note  also  that  the  parameter  estimates  cor¬ 
responding  to  the  coarsest  segmentation  differ  qualitatively  from  the  others.  This  segment 
length  is  so  long  that  intrinsic  variation  in  the  parameters  are  not  captured  appropriately. 
The  difference  between  the  parameter  estimates  for  the  finest  segmentations  is  mainly  due 
to  a  white  noise  estimation  residual.  For  the  log  intercept  process  this  residual  is  small 
and  we  have  essentially  obtained  what  we  sought,  a  parameter  estimate  that  is  stable  with 
respect  to  a  shortening  of  the  segmentation.  However,  for  the  slope  estimate  this  residual  is 
large.  This  leads  us  to  the  second  step  of  the  estimation.  In  this  step  we  carry  out  a  filtering 
procedure  in  order  to  remove  the  sampling  variability  that  is  segmentation  dependent  and 
particularly  strong  for  the  slope  process  p.  We  use  the  filtering  described  in  Section  4.4.3. 
From  the  first  step  we  obtained  the  estimates  (in  the  case  of  the  slopes) 

pl  =  pl  -f  wl  . 
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Figure  4.8:  Top  figure:  slopes  of  log-log  scale  spectra  from  wavelet  decompositions  based  on 
four  different  segment  lengths.  At  the  finest  resolution  the  scale  spectra  are  calculated  over 
nonoverlapping  segments  160m  long.  This  is  the  dotted  line  that  has  the  largest  variations. 
At  the  next  resolution  the  nonoverlapping  segments  are  327m  long  (plotted  with  a  dashed 
line).  The  other  two  resolutions  are  655m  and  1.31km  and  are  shown  with  dash-dotted 
and  solid  lines,  respectively.  The  bottom  figure  is  the  same  as  the  top  but  now  for  the  log 
intercept  of  the  scale  spectra.  The  gaps  correspond  to  inertial  ranges  with  fewer  than  five 
spectral  points. 


The  filter  computes  the  minimum  variance  unbiased  predictor  of  p1  given  p\ 

To  construct  the  minimum  variance  unbiased  filter  we  need  the  some  additional  estimates 
that  we  get  from  the  process  pl.  We  assume  that  the  slope  process  is  exponentially  correlated. 
We  need,  therefore,  its  correlation  length  lp  and  its  variance  ap 

E[(pi  ~  P)iPj  ~  P)}  =  exp(-L|i  -  j\/lp) 

with  L  the  length  of  the  segments.  Note  that  this  model  is  intrinsic  to  the  process  and 
does  not  depend  on  the  segmentation.  As  discussed  above  we  shall  model  the  sample  noise 
process  {m;}  as  white  and  we  only  need  its  level  cr-  .  We  estimate  the  parameters  ap,  lv  and 
using  the  empirical  variogram. 

For  a  time  series  X  =  (Xi)  of  size  N  the  empirical  variogram  with  lag  j  is  defined  by 

with  dependence  on  the  length  of  the  data  vector  N  not  shown.  Note  that  the  variogram 
is  essentially  unaffected  by  a  relatively  slow  variation  in  the  mean  P  of  the  process.  The 
mean  of  the  empirical  variogram  for  the  slope  process  P  is 

E[V{j)\  =  a2p(  1  -  exp(-L\j\/lp)  +  a2w 

and  from  this  we  obtain  the  parameter  estimates  by  fitting  to  the  empirical  variogram.  In 
particular,  afv  is  given  by  the  vertical  intercept  since  the  process  {u;*}  is  white. 

The  estimated  parameters  are  as  follows: 
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•  For  the  slope:  Vertical  intercept  cr^;  =  0.07,  horizontal  asymptote  level  cf2w  +  a2  — 
0.16,  and  correlation  length  800m  (five  segment  lengths). 

®  For  the  log  intercept:  Vertical  intercept  a2w  —  0.1,  horizontal  asymptote  level  a2,  + 
a2  —  3.0,  and  correlation  length  800m  (five  segment  lengths)  are  the  corresponding 
parameters. 

These  estimates  are  obtained  using  the  finest  segmentation,  160m,  and  only  the  data  in  the 
last  half  section,  the  high  turbulence  section.  The  estimates  were  based  on  slope  and  log 
intercept  processes  obtained  by  regression  as  above,  but  for  the  fixed  set  of  scales  7  to  12.  In 
this  case  the  magnitude  of  the  estimation  error  does  not  depent  on  the  segment.  We  show 
the  correlation  structure  in  Figure  4.9.  The  solid  line  is  the  fitted  model.  There  are  two 
parts  to  it,  the  sampling  noise  part  (the  intercept)  and  the  exponential  part  corresponding 
to  intrinsic  variation  in  the  parameters.  The  intercept  part  is  seen  to  match  well  with  that 
obtained  from  the  asymptotic  theory  presented  in  Section  6  which  predicts  its  value  to  be 
0.07  for  the  slope  and  0.05  for  the  intercept. 


Figure  4.9:  Variograms  for  the  fluctuations  of  the  slope  (top)  and  log  intercept  (bottom) 
processes  obtained  from  the  temperature  data.  The  stars  are  the  variograms.  The  solid 
lines  are  fitted  exponentials.  We  use  slope  and  log  intercept  estimates  from  the  data  with 
the  finest,  160m,  resolution. 

Using  the  estimated  model  parameters  we  carry  out  the  filtering  to  get  a  better  estimate 
for  the  slope  and  log  intercept  processes.  In  Figure  4.10  we  show  these  filtered  slope  and 
intercept  processes.  Note  that  after  the  filtering  all  three  segmentations  give  essentially  the 
same  result;  as  expected.  We  do  not  include  the  coarsest  segmentation  since  its  segment 
length  is  long  relative  to  the  interval  of  stationarity  and  intrinsic  variations  in  the  parameters 
are  being  smoothed  out. 

We  next  simulate  synthetic  temperature  profiles  using  the  parameters  estimated  above. 
Recall  that  we  only  aim  to  model  the  process  on  the  scales  corresponding  to  the  turbulent 
inertial  range.  We  therefore  scale  filter  both  the  aerotherman  data  of  Figure  4.1  and  the 
synthetic  temperature  profile  so  that  only  the  coefficients  in  the  Haar  wavelet  basis  corre¬ 
sponding  to  the  inertial  range  contribute.  Figure  4.11  shows  the  scale  filtered  processes. 
The  similarity  between  the  processes  is  striking.  Note  the  short  section  in  the  bottom  plot 
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Figure  4.10:  Filtered  slope  (top)  and  log  intercept  (bottom)  processes  of  the  actual  data  for 
the  three  finest  segmentations.  Note  that  the  filtered  processes  for  all  three  segmentations 
are  essentially  the  same.  The  dotted  line  corresponds  to  the  two  finest  resolutions.  The 
dash-dotted  and  dashed  curves  correspond  to  the  coarsest  segmentations.  The  filtering 
has  eliminated  the  differences  in  the  slope  and  log  intercept  processes  that  are  due  to  the 
segmentation  (as  in  Figure  4.8)  by  removing  the  white  noise  component  of  the  slope  and  log 
intercept  fluctuations  that  is  due  to  sampling.  Note  the  closeness  of  the  slope  estimate  to 
the  Kolmogorov  value  1.66  in  the  energetic,  second  half  of  the  data.  Note  also  the  increased 
value  of  the  log  intercept,  and  its  fluctuations,  in  the  energetic  second  half.  This  is  what  is 
expected  from  physical  considerations  (Washburn  et  al.  [34]). 


with  no  data  shown.  This  interval  corresponds  to  a  low  intensity  non-turbulent  interval 
in  the  atmosphere  where  a  laser  beam  propagates  without  much  distortion  due  to  medium 
heterogeneities.  Next,  we  assess  the  accuracy  of  the  parameter  estimates  by  applying  the 
estimation  procedure  to  simulated  data.  We  found  that  the  standard  deviation  between  the 
specified  and  estimated  parameters  for  the  local  slope  and  log-intercept  estimates  to  be  .15 
and  .3  respectively. 

4.6  Summary  and  conclusions 

We  have  analyzed  in  detail  estimation  of  fractional  Brownian  motion  based  on  the  scale 
spectrum.  Moreover,  we  have  generalized  the  estimation  procedure  to  a  local  power  law 
process.  For  such  processes  the  power  law  itself  (the  exponent  or  slope)  and  the  multiplica¬ 
tive  constant  (log  intercept  of  the  scale  spectrum)  are  not  constants  but  vary  slowly.  We 
estimate  the  slope  and  log  intercept  of  the  scale  spectrum  by  appropriately  segmenting  the 
data  and  then  removing  segmentation  effects  by  a  filtering  procedure.  The  slope  and  log 
intercepts  themselves  are  modeled  as  stochastic  processes.  The  estimation  of  the  power  law 
includes  the  identification  of  a  location  dependent  inertial  range,  that  is,  the  scale  range 
where  the  process  can  be  modeled  as  a  power  law. 

We  applied  the  estimation  procedure  to  an  important  set  of  aerothermal  data.  We 
found  that  the  average  of  the  estimated  slope  process  is  close  to  the  value  predicted  by  the 
Kolmogorov  scaling  theory.  To  get  a  faithful  representation  of  the  heterogeneity  in  the  data 
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Scale  filtering  of  data 


Figure  4.11:  The  actual  temperature  data  set  (top  figure)  obtained  by  synthesizing  it  from 
the  detail  coefficients  in  the  range  6  to  11  of  its  Haar  decomposition,  corresponding  to  about 
1  to  40  meters.  Information  from  this  range  of  scales  enters  into  the  estimates  of  the  slope 
and  log  intercept  processes.  The  bottom  figure  is  the  same  as  the  top  but  for  the  simulated 
data.  Note  the  striking  similarity  between  the  two  figures. 


it  is  however  important  to  incorporate  the  local  fluctuations  around  this  value. 

We  generate  synthetic  temperatures  in  two  different  ways.  First,  the  local  power  law 
process  is  simulated  conditioned  on  the  particular  parameters  estimated  from  the  data  (slope 
and  log  intercept).  Then  we  get  a  faithful  replication  of  the  (scale  filtered)  temperature  data. 
Second,  the  parameters  are  taken  as  random,  that  is,  sampled  from  their  model.  Then  the 
simulation  is  sample-wise  very  different  from  the  actual  data,  but  the  statistics  are  faithfully 
reproduced.  Such  realizations  can  be  used  to  define  synthetic  media  for  wave  propagation 
codes.  This  is  the  application  that  motivated  our  investigation.  A  complete  study  of  the 
effects  of  the  heterogeneity  in  local  power  law  model  on  wave  propagation  has  yet  to  be 
carried  out. 

An  important  aspect  of  the  model  that  we  use  is  separation  of  scales  in  the  variation  of 
the  estimated  parameters  (slope  and  log  intercept)  from  the  scale  at  which  we  sample  the 
process. 


4.7  On  estimation  of  fractional  Brownian  motion 


4.7.1  Statistics  of  wavelet  coefficients  for  fBm 

The  scale  spectrum  is  based  on  the  wavelet  coefficients.  We  now  consider  the  statistics  of 
the  Haar  wavelet  coefficients  for  fBm.  Next  we  will  analyze  the  scale  spectrum  and  power 
law  parameter  estimates  based  on  these  statistics.  We  assume  that  the  data  given  to  us  is 
the  wavelet  approximation  coefficients  at  level  zero.  That  is,  we  assume  that 


a0(n)  =  /  Bfj(x)dx  . 

Jn- 1 


(4.12) 
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The  wavelet  coefficients  are  therefore  normally  distributed  random  variables  with 


E[dj(n)\ 
Var[dj{n )] 

Cor[djl  ( ni)dj2 (n2)] 


Here 


1  =  2j2~j\  jl  <  j2 

(4.14) 

is  the  relative  scale, 

D  =  |(2n2-l)/-(2ni -1)| 

(4.15) 

is  the  relative  location  and 

6df{x)  =  f(x  +  d)~  2  f(x)  +  f(x  -  d) 

(4.16) 

=  0 

_  2  (l-2~2")  2j(2H+l) 

(2H  +  2)(2H  +  1) 

=  8(2^"-  i)l^/^|2g+2^Vi2/g|a!|2g+2>1^-  (4‘13) 


is  the  second  order  symmetric  difference.  These  results  are  derived  using  (4.1)  and  self¬ 
similarity  of  fBm.  Note  that  the  expression  for  the  correlation  of  the  wavelet  coefficients 
takes  on  a  universal  form  depending  only  on  the  relative  displacement  in  space  and  scale. 
The  derivation  of  the  precision  of  the  power  law  estimators  in  Section  4.7.3  is  based  on 
(4.13),  which  was  only  known  in  special  cases  before  (  Tewfik  and  Kim  [33];  Flandrin  [7]). 

Consider  the  correlation  of  wavelet-coefficients,  when 

j  =  \2i2(2n2-l) -2j'(2m -l)\2~h  +oo 


then 


tf/D  6i/D 
0 I/D )2  (1  /Dy- 


is  a  forth  order  central  difference  operator  and  its  expansion  gives  the  asymptotic  result 


Cord'd];} 


{H  +  1)(2H  +  \)H{2H  -  1) 
2(22H  -  1) 


If  the  scales  j\  and  jo  are  given,  then  for  no  ->  oo  and  rii  fixed  we  have  that  «  o'n'jV 

whereas  for  m  — >  oo  and  no  fixed  we  have  ^  ~  The  correlation  within  a  single  scale 
is  in  particular 


Cor[dfdf] 


(H  +  l)(2H  +  l)H{2H  -l) 
2{22H  -  1) 


2-2  H 

+  0 


For  general  wavelets  the  following  asymptotic  result  is  presented  in  Tewfik  and  Kim 
(1992)  and  Flandrin  (1992) 


E[cllld];]  ~  0(|2J‘2n2-2jlm|2(//-/?)) 

for  \2j2no  —  2^ln\\  —¥  +oo 

with  R  being  the  number  of  vanishing  moments  for  the  wavelet  (R  =  1  for  Haar  wavelets). 
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4.7.2  Statistics  of  the  scale  spectrum 

Using  the  above  results  for  the  wavelet  coefficients  we  get  a  characterization  of  the  statistics 
of  the  scale  spectral  points.  The  mean  is 


where 


E[Sj]  =  E[dj(n)2]  =  a2h(H)2j{2H+1) 

,  sin4(7r2J-1/) 


oc 


L 


o-'l/l 


(-xv-'fy- 


2  jdf 


h(H)  = 


(1-2~2H) 

(2H  +  2){2H  +  1)  ' 


The  normalized  variance  is 


Var{Sj] 

£[S;P 


,Era,m(Qmp 

'  NfE[Sj ]2 


1  f  4  /Vi_1 
=  Nj 


(4.18) 


(4.19) 


(4.20) 


with  Nj  =  2M~3}  the  number  of  detail  coefficients  at  level  j  and  2M  the  total  length 
of  the  data.  Note  that  here  and  in  the  sequel  we  suppress  the  dependence  on  M.  The 
matrix  C3  =  (C3nm)  in  (4.20)  is  the  covariance  matrix  of  the  wavelet  coefficients  at  scale  j, 
determined  by  (4.13),  and  pn{k)  is  the  corresponding  correlation  coefficient,  which  depends 
only  on  k  =  \ni  -  n2|  and  H.  Denote 


9 (R)  =  |  ^  53  ^  -  k)pH(k)  ~  2 1 

then 

VarjSj]  ~  g(H) 

E[Sj)2  ~  Nj 


(4.21) 


(4.22) 


for  Nj  large.  The  function  g(H)  is  computed  numerically  and  is  shown  in  Figure  4.12.  It 
depends  only  weakly  on  the  value  of  H ,  especially  when  H  <  1/2  which  is  the  case  of  interest 
to  us.  The  asymptotic  behavior  (4.22)  holds  only  when  the  correlation  p2H{k)  decays  to  zero 
in  an  integrable  way.  This  means  that  the  Hurst  exponent  must  be  restricted  to  H  <  3/4. 
In  the  Kolmogorov  case  H  =  1/3,  the  correlation  pH  decays  like 


Figure  4.12:  Plot  of  g{H)  showing  how  the  relative  variance  of  the  scale  spectrum  depends 
on  H. 
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Similarly,  if  H  <3/4,  the  normalized  covariance  can  be  expressed,  for  fixed  scale  sepa¬ 
ration,  as 


=  CovlS^Sj,]  En,m(^2)2  (. 

Un  ~  £&,]£&,]  NhNhE[Sh]E[Sh]  K  ^ 

with  CfJ'1*^2  being  the  cross  covariance  between  the  wavelet  coefficients  at  scale  j\  and 
determined  by  (4.13).  The  scale  spectral  points  5/  are  not  uncorrelated  for  different  j,  as  is 
sometimes  assumed.  For  the  power  law  parameter  estimation  we  note  that  for  fixed  j\  and 
3  2 


n  =  CovlSji ,  Sj2]  1 
nn  ~  E[Sh]E[Sh]  ~  JN~N~ 

as  Njx  and  Nj2  go  to  infinity.  These  results  can  be  derived  using  (4.13). 


(4.24) 


4,7.3  Fluctuation  theory  for  the  scale  spectra 

In  the  least  squares  fit  of  the  power  law  we  will  need  the  statistical  properties  of  the  logarithm 
of  the  scale  spectra.  We  summarize  the  relevant  facts  in  this  section,  with  the  details  given 
in  Appendix  1. 

In  the  large  Nj  limit,  the  distribution  of  the  scale  spectral  points  is  normal 


Si 

V3 


Af(0 ,g{H))  as  Nj  -4  oc 


with  Sj  =  E[Sj]  =  a2h(H)2^2HJrl\  and  where  h  and  g  are  defined  by  (4.21)  and  (4.19), 
respectively,  and  is  the  normal  distribution  with  mean  /i  and  variance  s 2 .  The 

derivation  of  these  results  requires  that  H  <  1/2.  The  covariance  of  the  fluctuations  is 


Cov  f 


Vi 


‘s/l Tj'y/Ni- 


—  D  ji 


(4.25) 


with  Dji  given  by  (4.24). 

The  fluctuations  in  Sj  are  small  and  in  Appendix  1  we  show  that  this  leads  to  the 
asymptotic  estimate 

l°g2(-Sj)  *  log2(^)+  jfi.  (4.26) 

VNi  ln(2) 

In  addition  to  showing  that  vj  is  asymptotically  normal,  we  also  show  in  Appendix  1  that 
for  any  j i  <  j-2  fixed  the  random  variables  Vj,  j\  <  j  <  jo  are  asymptotically  jointly  normal, 
as  Nj  -4  oo. 


4.7.4  Estimators  for  the  power  law 

Write  (4.26)  as 


log.)  (S'j)  «  \og.,(a2h(H))+j{2H+l)+  J  — 

V*j  ln(2) 


=  c  +  jp  + 


x/iVj  ln(2) 


(4.27) 
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with  the  slope  p  and  log  intercept  c,  to  be  estimated  from  data,  and  vj /  >/Nj  a  fluctuation 
term  that  is  characterized  in  the  large  Nj  limit  by  the  central  limit  theorem,  as  discussed 
in  the  previous  section.  In  view  of  the  above  we  can  then  use  generalized  least  squares  to 
estimate  b  =  [c,p]T  as 


b  =  ( XTD~lX)-lXTD-lY  (4.28) 

with  Y  —  [log2(SJ*1 ),  •  •  • ,  log2(5j2)]T,  and  where  j\ ,  •  •  ■ ,  jo  are  the  scales  in  the  range  under 
consideration, 

the  inertial  range.  The  dependence  on  j\ ,  jo  &  M  has  been  suppressed.  For  a  discussion  of 
generalized  least  squares  see  Hardle  [12].  The  design  matrix  X  is  defined  as 


X 


1  h 

1  jl  +  1 

I  h 


The  matrix  D  is  the  normalized  covariance  matrix  of  the  spectral  points  defined  in  (4.23). 
This  matrix  depends  on  the  value  of  H .  However  this  dependence  is  weak  so  we  can  use 
some  rough  estimate  of  H  that  is  usually  available,  like  H  =  1/3  for  the  aerothermal  data 
of  Section  5.  In  view  of  (4.27)  we  find  that  the  estimates  of  b  have  means 

E[c]  =  log2(c72Mtf))  (4.29) 

E\p]  =  277  +  1 


and  covariance 


Cb  =  (XtD~1X  ln(2)2)-1 

ln(2)"2 


E  E  DTjl  -ZZDji 
-EE D$i  ZZDjij 


.  [E E E E DJj-  * i-(EE Dl?)2] 

for  large  Njx .  From  the  above  we  arrive  at  the  sought  after  parameter  estimates 

H  =  (p  — 1)/2 

log2(<J2)  =  c-log2(M^))- 


(4.30) 

1 

Xu 


(4.31) 

(4.32) 


We  next  analyze  the  precision  of  these  estimators.  We  find  that  the  estimates  H  and 
log2(<r2)  are  normally  distributed  with  variances 

Var[H }  »  Cb(2,2)/4  (4.33) 

Har[log^r2)]  «  C'b(l,  1)  —  2Cb(l,2)  log2(h(H))'  +  Cb(2, 2)  [\og2(h(H))f  (4.34) 


for  Njx  large  with  Cb  defined  in  (4.30)  and  h  in  (4.19).  The  variance  of  the  estimators  is  of 
order  1  /Njx. 

Long-memory  processes  can  also  be  modeled  via  state  space  models.  Maximum  likelihood 
estimators  (MLE)  for  such  models  are  analyzed  in  Dahlhaus  [4]  and  Fox  and  Taqqu  [8]. 


4.7.5  Illustration  of  precision 

The  above  result  on  the  distribution  of  the  estimators  may  be  validated  by  numerical  sim¬ 
ulation.  We  generate  synthetic  realizations  of  fractional  Brownian  motion  with  known  pa¬ 
rameters.  Then  we  use  the  above  algorithm  to  estimate  these  parameters  and  compare 
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the  precision  of  the  estimates  with  the  predicted  precision.  Actually,  we  do  not  simulate 
fractional  Brownian  motion  but  rather  the  observed  sequence  (4.12),  the  approximation  co¬ 
efficients  at  level  zero.  The  synthetic  realizations  are  generated  by  a  simulation  algorithm 
similar  to  the  one  discussed  in  Omre.  S0lna  and  Tjelmeland  [23]. 

In  Figure  4.13  we  compare  the  asymptotic  law  of  the  estimators  with  simulations.  It 
is  clear  that  the  theoretical  normal  distribution  predicts  accurately  the  distribution  of  the 
estimators. 


Log  intercept  estimates  Exponent  estimates 


log  intercept  exponent 


Figure  4.13:  Distribution  of  power  law  parameter  estimates.  The  parameters  were  estimated 
based  on  realizations  of  length  210  with  H  =  1/3  and  a 2  =  2.  The  power  law  fitting  is 
done  over  scales  1---5.  The  solid  line  is  the  empirical  distribution  associated  with  1000 
realizations,  the  dashed  line  is  the  theoretical  distribution  and  the  vertical  lines  are  the 
specified  parameters. 

We  show  how  the  variance  of  the  slope  estimate  p  depends  on  the  number  of  scales  in 
the  linear  fit  and  the  total  number  of  data-points  in  Figure  4.14.  The  dotted  and  solid  lines 
correspond  to  using  respectively  scales  1,2,3  and  scales  1,2, 3, 4, 5  and  Hurst  parameter 
H  —  1/3.  The  crosses  give  the  Cramer- Rao  lower  bound  (Abry  et  al.  1995;  Wornell  and 
Oppenheim  1992):  1/(2M+2  ln2(2))  for  the  case  that  the  wavelet  coefficients  are  uncorre¬ 
lated.  It  is  seen  that  only  a  few  scales  are  needed  for  the  variance  of  the  slope  estimator  to 
be  close  to  the  bound. 


Figure  4.14:  The  variance  of  the  slope  estimate,  computed  by  (4.33)  with  H  =  1/3,  is 
plotted  as  a  function  of  the  number  of  data  points  2M .  The  dotted  line  is  for  an  estimate 
with  three  scales,  the  solid  line  with  five  scales.  The  crosses  correspond  to  the  Cramer- Rao 
lower  bound.  Note  that  even  when  the  number  of  data  points  is  not  very  large  the  estimate 
of  the  slope  does  not  depend  sensitively  on  the  number  of  scales  used. 

Note  that  our  objective  is  estimation  of  H  when  H  as  1/3  and  that  our  asymptotic 
result  is  valid  only  when  H  <  1/2.  Several  authors  (Ninness  1998)  have  observed  that  slope 
estimators  degrade  when  H  become  large  ( H  ss  1).  This  is  because  in  that  case  the  wavelet 
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coefficients  exhibit  long  range  correlations. 

In  the  application  of  the  estimation  of  the  power  law  parameters  to  the  aerothermal  data 
in  Section  5,  which  we  model  by  a  local  fractional  Brownian  motion,  we  have  to  account  for 
spurious  noise  effects.  This  noise  affects  the  small  scales  mostly.  We  can  make  the  power 
law  parameter  estimates  more  robust  by  letting  the  weight  matrix  D  in  (4.28)  depend  on 
an  additive  noise  parameter  which  must  be  adjusted  appropriately  from  rough  prior  data 
analysis.  This  is  done  in  Section  4.5.2  without  a  detailed  discussion. 


4.8  Central  limit  theorem  for  scale  spectra 


In  this  appendix  we  derive  a  central  limit  theorem  for  the  scale  spectral  points,  Sj,  and  also 
the  log  transformed  spectral  points  log2  (Sj),  where 


Sj  = 


Nj 


with  the  Haar  wavelet  coefficients  of  fractional  Brownian  motion  at  scale  j  with 

H  <  1/2.  The  total  number  of  data  points  is  2M  and  Nj  =  2M~j .  We  will  consider  the 
limit  Nj  — ¥  oo.  This  means  that  M  is  large  and  that  j  is  not  too  close  to  M,  that  is,  the 
central  limit  theorem  will  not  be  valid  for  the  large-scale  spectral  points  Sj  with  j  %  M. 

We  will  use  the  following  version  of  the  Berry-Esseen  theorem  (Feller  [6]). 

Let  the  yi  be  independent  random  variables  such  that 

EM  =  o,  E[yf\  =  al  E[\y*\ }  =  Pi 


and  define 


Sn  =  rn  =  X>- 

z=l  i—  1 

Let  Fn  be  the  distribution  of  the  normalized  sum  yi/sn. 
Then  for  all  x  and  n 


\En{x)  —  N(x)\  <  6rf  (4.35) 

where  N  stands  for  the  mean  zero,  unit  variance  normal  distribution. 


4.8.1  Central  limit  theorem  for  Sj 

To  use  this  theorem  we  will  first  transform  the  sum  of  squares  of  wavelet  coefficients  to  a 
sum  of  independent  random  variables.  Denote  the  vector  of  wavelet  coefficients  at  scale  j 
by  $  =  [dj(  1),  •  •  • ,  dj(Nj)]T  and  the  covariance  matrix  associated  with  cP  by  CF  Then  d i 
has  the  same  law  as 

and  in  law 


Nj 


Ni 


Sj  _  "  ST+ AlASAi(,?? 


1)] 


(4.36) 
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where  iji  are  independent  zero  mean  and  unit  variance  normal  random  variables,  Sj  =  E[Sj] 
and  A:-  are  the  eigenvalues  of  Cj /Sj,  that  is  the  correlation  matrix  of  dj . 

Let 


Yj  = 


1  N’ 

75  & 


with  yj  —  Aj (rjf  -  1).  Note  that  the  yj  are  independent  random  variables  such  that 
E[y{]  =  0,  E[(y{f]  =  2(Aj)2,  E[\(y{)3\}  <  28(A3)3  . 

According  to  (4.35)  we  need  only  show  that 


J  s  eSaMI 

is  small  for  Nj  large. 

Using  (4.20)  we  find  that 

AO 

EM"  =  H^nm/Sj)2  =  N2 V ar[Sj }/(2S2)  »  Nj  g(H)/2  . 

i—  1  n,m 

Below  we  show  that 

M  <  K 

for  some  constant  K  >  1  independent  of  Nj .  Hence 


(4.37) 


(4.38) 


(4.39) 


E(A])3  <  NjK3  .  (4.40) 

1=1 


From  (4.38)  and  (4.40)  we  find  that 


jj 


(2  /r2/g(g))3/2 


and  goes  to  zero  as  Nj  oo.  From  the  Berry-Esseen  theorem  we  conclude  that  the  distri¬ 
bution  of 


Yj/y/fKH) 


E 


tends  to  the  standard  normal  distribution  as  Nj  becomes  large. 
Thus,  for  Nj  large 


Sj 


Sj{  1  +  Cj 


in  law,  with  tj  a  standard  normal  random  variable. 

Let  us  now  show  that  the  estimate  (4.39)  is  valid.  The  diagonal  entries  of  the  correlation 
matrix  are  all  equal  to  one.  By  the  Gersgorin  disc  theorem,  |A^  —  1|  <  \^in\/Cn.  But 

the  sum  remains  finite  as  Nj  oo  because  we  have  assumed  that  H  <  1/2  and  so  we  can 
use  (4.17).  Rewriting  this  inequality  with  Aj  on  the  left  side  we  get  (4.39). 
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4.8.2  Central  limit  theorem  for  log2(Sj) 

We  consider  next  a  central  limit  theorem  for  log2  [Sj].  As  in  (4.36)  we  write  the  scale 
spectrum  in  terms  of  centered  random  variables  Yj 


Sj  —  *Sj(l  + 


JL) 


Taking  logs  and  expanding  in  a  Taylor  expansion  with  remainder  we  get 


log-2  (Sj) 


logaC^j)  + 


*7 


2£(Yj)2Nj  In  2 


with 


<  1,  1  +  Yj/y/N]  > 

<i  +  yi/y^7,i> 


when  Yj  >  0, 
when  Yj  <  0. 


We  showed  above  that  the  central  limit  theorem  holds  for  Yj.  We  show  below  that  Y2 / (2£(Y))2) 
is  0(1)  as  Nj  oo.  We  therefore,  have  a  central  limit  theorem  for  log2(Sj)  as  well. 

By  the  Cauchy-Schwarz  inequality 


E 


Y2 

I — 1—1 

k(ij)2 


We  show  that  Jj  =  0(1)  as  Nj  oo.  Consider  first  E[Y )4].  Using  (4.39)  we  find 
E[Yf]  =  g[(E^^~1)])4] 

<  K*(NjE[(rfi  -  l)4]  +  Nj  (Nj  -  l)6E?[(rft  -  1  )2})/N]  =  a  +  c2/Nj 
with  C{  constants.  Next  consider 


Ij  =  EiaYj)-*]  <  E[aYj)-4I{Yi<0)}  + 1 

<  £[(i  +  w^)-Vi<o,]  +  i 

with  I A  the  indicator  function  of  the  set  A.  In  order  to  bound  Ij  we  replace  1  +  Yj/  y/N] 
by  Zj  such  that  w.p.l  Zj  <  (1  +  Yj/y/Nj). 

Note  that 


i  +  W*7  =  ^4 

and  that  A]  =  Nj,  which  is  the  trace  of  the  correlation  matrix  for  dY  Let  q  be  the 
number  of  eigenvalues  that  exceed  1/2.  In  view  of  (4.39)  we  find 

q  >  Nj  =  Nj(2K  —  1)_1  =  NjK  . 

Thus,  we  can  define  Zj  as 


1E&  m  =  !EM 

2  Nj  2  Nj 
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Hence 


Ij  <  E[Z~4I{Y.<0)]  +  1  <  (2/ K)~4  N4  E[Z~4]  +  1 

with 


2-1 


in  law.  The  random  variable  Zj  has  law  given  by  the  Gamma  density 

7(cmo(z)  =  exp(-Qx) 

with  parameters  a  —  1/2  and  i/  =  iVj/2.  It  follows  that  for  Nj  >  10 


E\Z-4)  =  (<*rrfr-4)  =  _ I _ 

1  '  (q)— 4r(j/)  (Nj  -  2)(Nj  -  A)(Nj  -  6)(Nj  -  8) 

Hence,  Ij  <  Kn  and 

Jj  <  y/ K" (a /Nj  +  a)  <  K 


for  Nj  >  10/ K. 

From  the  above  we  conclude  that  for  Nj  large 


log  2(Sj)  «  log  2(Sj)+€j 
in  law,  with  tj  a  standard  normal  random  variable. 


9(H) 
Nj  ln(2) 


4.8.3  Joint  density  of  spectral  points 

We  show  that  the  spectral  points  are  asymptotically,  as  Nj  — >  oo,  jointly  normal.  From 
this  follows  that  the  central  limit  theorem  holds  also  for  the  power  law  parameter  estimates. 
Consider  the  distribution  of  Y  =  axSjx  +  a2Sj2  with  ax  ^  0  k  a2  ^  0  and  j\  <  j2  fixed 
parameters.  The  central  limit  theorem  for  this  quantity  can  be  shown  essentially  as  in 
Section  4.8.1.  Let  d-71  and  dj-  be  the  vectors  of  wavelet  coefficients  at  the  two  scales.  Define 


and  d  =  Dd  with  D  a  diagonal  matrix  whose  first  Njl  elements  equal  ax  and  last  Nj. 
elements  equal  a2  2 J2  “-71 .  Then 


Y  =  drd/NJX  . 


By  a  transformation  to  independent  variables  and  an  application  of  the  Berry-Esseen  theo¬ 
rem  we  find,  as  in  (4.37),  that  we  need  to  show 


J 


Nji+Nfr 


(AO3 


Ei=i+iVj2(^i)2]3/2 


(4.41) 
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is  small  for  Njx  large.  Here  A  *  are  the  eigenvalues  of  Cov(d)  =  DCD  with  C  the  covariance 
matrix  of  d.  Note  first  that  from  (4.20)  it  follows  that 

+Nj2 

£  A*  >  (ai)2E(^m)2  «  («i)2^-,(4)2^)/2. 

i—  1  n,m 

Hence,  (4.41)  follows  if  we  can  bound  A ;  uniformly  with  a  bound  independent  of  Njv .  We 
find,  using  (4.13)  and  the  Gersgorin  disc  theorem,  that  for  some  constant  L  and  1  <  j  <  Nj: 

|Aj  ~  ai§ji\  ^  LSji 

since  the  rows  of  C  are  absolutely  and  uniformly  summable  with  a  bound  that  is  independent 
of  Njl .  Hence  A;  <  LSjl  for  some  constant  L.  A  similar  argument  holds  for  Njl  <  j  < 
Njl  +  Nj2  .  Therefore,  (4.41)  and  the  central  limit  theorem  for  Y  follows. 
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Chapter  5 

Spectral  Analysis 


5.1  Introduction 


For  a  segmentation  with  segment  length  T  and  sample  interval  A t,  the  measurable  frequency 
range  is  2ir/T  <  u  <  ir/At,  which  corresponds  to  the  scale  range  T  to  A £.  While  A  is  usually 
dictated  by  the  intrinsic  bandwidth  of  the  process  being  sampled  or  the  instrumentation, 
the  choice  of  T  has  always  been  somewhat  arbitrary.  The  SpecLab  procedure  relies  on 
matching  scale  spectra  to  a  power-law  model.  In  principle,  standard  spectral  analysis  could 
have  been  used,  but  there  are  compelling  reasons  to  use  scale  spectra  for  this  purpose.  The 
material  in  this  chapter  has  two  objectives.  The  first  is  to  demonstrate  the  problems  that 
are  encountered  when  standard  spectral  analysis  procedures  are  applied  to  fBm.  The  second 
is  to  provide  guidlines  for  translating  the  SpecLab  structure  function  and  Hurst  exponent 
to  spectral  strength  and  power-law  index. 

To  introduce  the  analysis  consider  a  stationary  power-law  process  that  has  an  outer  scale. 
In  that  case  both  the  autocorrelation  function  and  the  spectral  density  function  are  well 
defined.  Furthermore,  the  process  itself  admits  a  spectral  decomposition.  There  is  a  wealth 
of  material  that  describes  standard  analysis  procedures  that  can  be  applied  using  discrete 
Fourier  transforms  (DFTs)  to  estimate  the  spectral  density  function  (SDF).  For  reference, 
the  essential  elements  of  this  material  for  our  purposes  are  summarized  in  Section  5.2.  The 
development  in  Section  5.2.2  shows  that  the  structure  function  of  a  stationary  power-law 
process  admits  the  following  limiting  form  for  large  outer  scale: 


D{y) 


C  r((3-^)/2)22~‘- 
V?  (t>-i)r(i//2) 


=  o'2  |j/| 


V- 1 


(5.1) 


The  relation  is  valid  for  1  <  v  <  3  over  the  y  range  of  scales  that  encompass  the  power-law 
segment.  The  unbounded  power-law  form  of  D(y)  does  not  have  a  Fourier  transform,  but 
(5.1)  implies  the  following  scale-free  relation  between  the  structure  constant  a2  and  the 
spectral  strength  parameter  C : 


C  F((3  —  v)/2)22~ 

~  \fz  (v  -  i)r(i//2)  • 


If  the  process  is  obtained  from  a  probe  sampling  a  two-  or  three-dimensional  process 
via  translation  and  path  integration,  additional  scale-free  transformations  relate  the  one¬ 
dimensional  parameters  to  an  appropriate  scale  free  characterizations  of  the  actual  medium. 

In  the  remainder  of  this  Chapter  we  first  calculate  the  bias  characteristics  of  the  standard 
periodogram  estimator  when  applied  to  fBm  as  opposed  to  a  realizations  of  a  stationary 
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process.  The  results  show  that  the  non-stationary  characteristics  of  fBm  introduce  a  more 
severe  bias  that  is  encountered  with  stationary  processes.  Furthermore,  that  bias  will  defeat 
attempts  to  estimate  the  spectral  slope  independent  of  the  spectral  strength.  We  then  show 
that  for  the  Haar  wavelet  applied  to  fBm  the  relation  (5.2)  is  exact.  The  combination  of 
these  two  observations  strongly  support  the  wavelet-based  analysis. 


5.1.1  Periodogram  Bias  for  fBm 

Consider  the  periodogram  estimator  Pn  —  \zn\ 2  / N ,  where  zn  is  the  discrete  Fourier  trans¬ 
form  of  a  stationary  process.  The  development  in  Section  5.2.2,  suggests  that  E  [Pn]  « 
Cu~u / At.  To  test  this  conjecture,  the  expectation  of  the  periodogrm  can  be  computed 
from  the  defining  properties  of  fBm.  With  some  straightforward  manipulations,  it  follows 
from  the  definition  of  the  DFT  and  the  fBm  covariance  function  that 


E 


\z(n)\-/N 


=  (e£Lo‘ 

4 -a2  (1  —  lml  /N)  m'2H  cos(2mnm/N). 


(5.3) 


Figures  5.1  shows  the  result  of  averaging  10  periodogram  estimates  from  independent  2 12 
point  fBm  realizations  with  H  —  1/3.  Superimposed  on  the  averaged  periodograms  is  the 
expectation  derived  from  (5.3)  with  a2  —  2~2dH ,  where  d  is  the  dyadic  dimension  and  Sn  is 
the  Dirac  delta  function.  The  nominal  power-law  with  C  derived  from  (5.2)  has  also  been 
overlaid  on  the  periodogram  plot.  We  note  first  that  the  averaged  periodogram  agrees  well 
with  the  expectation  up  to  the  highest  frequencies.  The  departure  comes  from  a  small 
modification  of  the  fBm  process  introduced  to  represent  integration  over  the  sample  period 
rather  than  instantaneous  sampling.  More  important,  although  the  expectation  appears  to 
be  reasonably  power-law  in  form,  the  deviation  from  ideal  power-law  is  frequency  dependent 
and  varies  with  the  Hurst  exponent.  This  behavior  would  not  be  expected  for  stationary 
processes  as  summarized  in  Section  in  5.2.2. 


5.1.2  Wavelets  and  Scale  Spectra 


Wavelet  analysis  has  received  considerable  attention  over  the  past  decade,  particularly  with 
the  discovery  by  Mallat  [38]  of  the  intimate  relation  between  multiresolution  filtering  and 
wavelet  decompositions.  The  book  “Wavelets  and  Filter  Banks”  by  Strang  and  Nguyen[39] 
covers  the  material  in  detail. 

Continuous  Wavelet  Transforms 

The  continuous  wavelet  transform  (CWT)  of  x(t')  is  defined  as 

r°°  1  tf  -  r 

Fw(s,t)  =  /  x{t')—=w{ - )eft\  (5.4) 

J-OO  V*  S 

where  w(t)  is  the  “mother”  wavelet.  The  parameter  s  is  a  scale  parameter  akin  to  period 
(1//)  while  r  locates  the  region  where  the  scale  information  is  extracted.  To  interpret 
(5.4),  assume  that  x(t)  admits  a  spectral  representation  of  the  form 

x(t)  =  J  \JQ{uj)  exp{iut}d£(uj),  (5.5) 
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Hurst  Exponent  =0.33333 


Figure  5.1:  Average  of  periodograms  from  10  independent  realizations  of  fBm  for  H  =  1/3. 
The  magenta  curve  is  the  expectation  value  of  each  frequency.  The  blue  curve  is  the  ideal 
power-law. 


where  d((u)  has  the  formal  “white-noise”  property  E  [d((uj)dC  (a/)]  =  8(uj  —  u/)/(27t).  The 
well-known  Fourier-transform  relation 

E[x(t)x*(t1)}  =  J  Q(uj)exp{iuj(t  -  (5.6) 


follows  from  the  white-noise  property.  Substituting  (5.5)  into  (5.4)  gives  an  equivalent 
spectral-domain  representation  of  the  CWT: 


/oo 

\f sQ(lj)w* (suj)  exp{icjT}d((ui). 

-OO 


(5.7) 


Exploiting  the  white-noise  property  as  before  generates  the  r-independent  scale  spectrum 


E  |iv,(s;r)|2]  =  J  Q(u)s\w(su)\2^. 


(5.8) 


For  a  power-law  process  spectrum 


E 


|w>MI 


•2  du 
2tt 


Csv. 


(5.9) 


Since  s  oc  1/cj,  the  wavelet  scale  spectrum,  unlike  the  fft-based  estimator,  evidently  extracts 
the  correct  power-law  form.  Wavelets  have  the  property  that  they  integrate  to  zero,  which 
implies  that  lim  w(scj)  =  0.  Thus,  with  an  appropriate  wavelet  choice,  the  bias  term  may 

oj — >-0 

remain  finite  in  spite  of  the  singular  behavior  of  the  power-law  u  =  0.  Ninness  [40]  and  Abry 
et  al.[l]  have  shown  independently  that  the  expectation  of  the  scale  spectrum  for  an  fBm 
process  does  have  the  scale-free  property  implied  by  (5.9). 
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Discrete  Wavelet  Transforms 

The  discrete  wavelet  transform  (DWT)  is  obtained  by  taking  t  =  kAt  and  s  =  2^”1A £, 
whereby 

i  r°° 

dJk  =Fw{V^At,kAt)  =  . c.  - /  x{t*)w((t9 / At  -  k)/2^1)  dt'.  (5.10) 

v  2J  i  At  J— oo 

It  is  very  inefficient  to  evaluate  (5.4)  by  direct  numerical  integration,  but  following  Mal- 
lat’s  [38]  discovery  of  the  intimate  relation  between  multiresolution  filtering  and  wavelet 
decompositions,  the  DWT  can  be  computed  efficiently  without  explicit  use  of  the  wavelet 
or  the  scale-function  that  generates  it.  Two  sets  of  filter  coefficients  are  applied  recursively 
with  downsampling  to  generate  a  wavelet  scale  decomposition  that  approximates  (5.4)  at 
logarithmically  spaced  scales.  The  results  summarized  here  are  developed  in  detail  in  the 
textbook  by  Strang  and  Nguyen  [39]. 

Each  wavelet  has  characteristics  that  strongly  affect  its  applications.  The  MatLab 
Wavelet  Tool  Box  provides  a  family  of  wavelets  that  are  defined  by  orthogonal  low-  and  high- 
pass  filter  coefficients  Lo.Dm  and  HiJDm  that  can  be  generated  by  the  function  orthfilt . 
Both  filters  have  L  wavelet-dependent  coefficients.  The  downsampled  low-pass  filter  output 
at  level  (scale)  j  generates  approximation  coefficients  ak  while  the  high-pass  filter  generates 
the  corresponding  detail  coefficients  dk  that  form  the  DWT.  The  recursion  is  initiated  at 
level  0  with  the  data  samples: 

a°k  =  xk  for  k  =  0, 1,  •  •  • ,  N  -  1,  (5-11) 

With  N  =  2m  for  j  =  1,  •  •  • ,  M,  the  recursive  filtering  and  downsampling  operations  are 
defined  as 


ak  “*  S  L°~D rn ^q(m)  (5.12) 

771  =  0 

4  =  E  (5.13) 

771  =  0 

where  q(m)  =  2k  -  (m  -  1)  +  a  identifies  the  scale  sample  from  the  previous  level  that  is 
multiplied  by  the  mth  filter  coefficient.  The  2k  increment  accounts  for  the  downsampling. 
The  offset  a  is  a  phase  factor  that  determines  where  the  current  filter  output  is  referenced. 
The  recursion  is  not  completely  specified  until  the  procedures  for  accommodating  the  first 
L  —  1  and  the  last  L  —  1  samples  are  specified.  The  default  MatLab  procedure  is  to  affix 
L  —  1  leading  and  trailing  zeros  to  the  data  set.  As  a  consequence,  the  approximation  and 
detail  coefficient  vectors  at  level  j  —  1  are  of  length  Ld(  1)  =floor[(2M  +  (L  -  l))/2] ,  where 
floor[.]  means  the  largest  integer  not  exceeding  the  number  in  the  brackets.  The  length  of 
the  detail  coefficient  vector  for  j  >  1,  as  can  be  verified  upon  executing  the  MatLab  Wavelet 
Tool  Box  wavedec  routine,  is 

LD(j)  =  floor  [( 2m  +  (L-  1)(2j  -  l))/2>']  -  (5.14) 

Note  that  Lo(j)  >  LmU)  =  2M~J,  where  L,\f(j)  is  the  nominal  length  of  the  data  vector 
after  downsampling.  The  orthogonality  of  the  filters  implies  the  norm  preserving  relation 

N- 1  LD{j)  t>  j  Lo{rri) 

£  M  +I IS  (dk)2-  (5-15) 

k= 0  k=0  m=l  k= 0 
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The  wavelet  scale  spectrum  is  obtained  by  summing  the  intensity  of  a  truncated  set  of  the 
coefficients  at  each  scale: 

1  LdU)  /  \  2 

*  ■  £/*  («  ■  (516) 

with  UJk  symmetrically  zero  at  the  leading  and  trailing  edge  of  the  dJk  vector  to  reduce  its 
contributing  length  to  L.u(j)- 


5.1.3  The  Haar  Wavelet 


Because  the  properties  of  fBm  are  defined  in  terms  of  its  covariance  function,  we  proceeds 
from  relations  such  as  (5.10)  or  (5.12)  and  (5.13).  Closed-form  solutions  are  difficult  to 
come  by,  but  the  Haar  wavelet  has  a  particularly  simple  form: 


/,v  _  f  1  forO  <  t  <  1/2 

WHaar(t)  -  j  _j  <  t  <  1  - 

Substituting  the  Haar  wavelet  into  (5.4),  for  example,  yields 

1  fs/ 2 

F-w(s,t)  = -yz  /  Wu  +  r)  -  x(a  +  r  +  s/2)]  du. 

Vs  Jo 

Computation  of  the  expectation  of  the  scale  spectrum  is  straightforward: 


(5.17) 


(5.18) 


E  |F„(5,r)|5 


+  \u  —  u'  +  1/2\2H  —  2  | u  —  u,|2'W|  dudu'  s2H+l 


(5.19) 

J 

Carrying  out  the  integrations  and  using  the  equivalence  between  (5.10)  and  (5.13)  shows 
that 


E 


di 


1-2 


-2  H 


One  can  also  show  that 


(2H  +  2)(2H  +  1) 
sin4(w/4) 


2j(2H+l) 


i^Mr  = 


(u/4  )2 


(5.20) 


(5.21) 


Thus,  from  (5.9)  the  following  equivalent  expression  involving  the  spectral  strength  param¬ 
eter  is  obtained: 

,4/ 


E 


d{ 


1  =  P  Pjp-  ?)  .  (5.22) 

J  Vi-oo  (w/4)2  2ttJ 


Equating  (5.20)  and  (5.22)  gives  an  exact  relation  that  defines  a2 /C: 


C 


(£ 


w 


~(2W+i)  sin4(u;/4)  du\  _  2  1-2  ‘ 


(w/4)2  2ir 


1-2  H 


=  er 


(2H  +  2)(2H  +  1)' 


(5.23) 


Mathematica  (courtesy  of  Raul  Martinez  at  Vista)  produced  the  following  solution,  which 
is  evidently  valid  at  H  =  0: 


2(2~2H 


1/12 

1)  cos(if7r)r(-2(l  +  H)) 


H  =  1/2 
0  <=  H  <  1/2 
1/2  <  H  <  1 


(5.24) 
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Figure  5.2:  Average  of  wavelet  scale  spectra  from  10  independent  realizations  of  fBm  for 
H  =  1/3.  The  magenta  curve  is  the  expectation  value  for  each  scale  converted  to  frequency. 
The  blue  curve  is  the  ideal  power-law. 


Substituting  from  (5.2)  for  a2/C  provides  what  we  might  expect  to  be  an  approximate 
relation: 

I  .-(2/r+i)  sinV/4)  dw\  r(l-H)2~2H  1-2~2H 

1  1  (w/4 )2  2~)  y/^HT(H  + 1/2)  (2H  +  2)(2H  +  l)'  K  ’ 

Figure  5.2  shows  the  raw  scale  spectra  derived  from  fBm.  The  coefficients  are  plotted 
as  red  asterisks  at  the  frequencies  uJNyquist/2j~l  ■  The  solid  red  curve  shows  the  theoret¬ 
ical  expectation  value  derived  from  (5.20).  The  blue  asterisks  have  the  bias  removed  by 
using  (5.23).  The  continuous  blue  curve  is  the  theoretical  SDF.  In  sharp  contrast  to  the 
periodogram  estimator,  the  scale  spectrum  captures  the  power-law.  The  offset  from  the 
ideal  power  law  depends  only  on  the  Hurst  exponent  and  is  easily  compensate.  The  scale 
spectrum  removes  the  noise  fluctuation  without  biasing  the  power-law  index.  Note  that  the 
norm  preserving  relation  is  of  little  use  in  calibrating  the  strength  scale  because  it  contains 
a  non-stationary  contribution  from  the  average  wavelet  coefficients. 
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5.2  Background 


The  following  material  summarizes  basic  relations  that  are  used  in  ScaleSpec. 

5.2.1  Discrete  Fourier  Transforms 

For  a  sequence  of  complex  numbers  x(k)  where  0  <  k  <  N  —  1  the  discrete  Fourier  transform 
(DFT)  and  its  inverse,  as  evaluated  by  the  MatLab  functions  fft  and  ifft  are  defined  as 
follows: 


J(n)  —  Hk=o  x(^)  exp{— 2~ink/N}  (5.26) 

x{k)  =  jj  Zn=o  %{n)  exp{2~ink/N)  (5-27) 

The  Parseval  relation 

AT-i  N-i 

E  mi2  =  n  E  ijwi2  (5-28) 

k=0  n- 0 

is  readily  derived  from  (5.26)  or  (5.27). 

If  x(t)  and  x(uj)  are  related  via  the  continuous  Fourier  and  inverse  Fourier  transforma¬ 
tions 


x(uj)  =  x(t)  exp{—iut}dt  (5.29) 

x(t)  =  f™xx(oj)exp{iwt}^,  (5.30) 

it  can  be  shown  that  xp(kAt)/ At  and  xp(nAuj)(2ir  /  Au)  form  an  exact  discrete  Fourier 
transform  pair.  The  periodic  or  “aliased”  form  of  the  continuous  functions  are  defined  as 

xp{x)  =£m=-oo  x(t  +  mT)  (5.31) 

=  Em=— oo  +  2xm/At)  (5.32) 

where  At  =  T/N  and  Acu  =  2i r/T,  whereby  AtAoo  =  27r/Ar.  The  sampled  functions 
themselves  are  the  tapezoidal-sum  approximations  to  the  corresponding  integrals.  For  a 
continuous  Fourier  transform  pair,  one  chooses  At  and  N  so  that  the  intervals  —  T /2  <  t  < 
T/2  and  —n/At  <  u  <  tt/A t  contain  the  essential  content  of  the  corresponding  function  and 
its  Fourier  transform.  The  fact  that  natural  reference  for  the  functions  is  t  =  0  or  u  =  0, 
whereas  the  discrete  transforms  place  zero  at  the  beginning  of  the  interval  establishes  the 
need  for  a  cyclic  shift  routine,  which  MatLab  provides  as  fftshift. 


5.2.2  Stationary  Processes 

A  stationary  stochastic  process  admits  a  spectral  representation  of  the  form 


/oo 

v/CM  exp  {iwi}  dCM, 

-OO 


(5.33) 


where  d((u)  has  the  formal  “white-noise”  property  E[d((cj)  d^fa')}  —  S(co  —  o/)/(27t ). 
With  t  =  kAt ,  it  follows  from  (5.26)  that 


x{nAuj) 


1  —  exp  {iT(uj  -  nAu)} 
1  —  exp{iAf(o;  —  nAui)} 


rfC(w). 


(5.34) 


51 


The  raw  periodogram  estimator  is  defined  as 

Pn  =  \in\2/N.  (5.35) 

The  expectation  of  the  periodogram  is  readily  computed  by  using  the  white-noise  property: 

E 

L'  "■  *  j  ,/~uu 

w  Q(nAu)/ At.  (5.36) 

The  approximation  holds  for  large  N.  Note  that  DFTs  satisfy  the  Parseval  relation 


'  |2  /,vl  _  _L  r°°  D(,  ,’1  sin2(f  (u-«Am))  du 

Cn  1  'l  \  -ivJ-ooWsin2^(u_nia))2iI 


1  N—l  N- 1 

j\7  £  i^i2=  E  w2- 


N  n 

n-0  k— 0 

Upon  taking  expectations  and  using,  (5.36)  we  find  that 


(5.37) 


N- 1 

Q(n&u)/At  =  N  (x2) 

n=0 


(5.38) 


or 

E  Q(nAu)^r  =  (x2)  •  (5.39) 

n— 0 

That  is,  in  expectation  for  large  N  the  sum  of  the  periodogram  estimator  is  the  discrete 
approximation  to  the  corresponding  integral  relation  for  the  SDF  of  a  stationary  process. 

Power  SDFs  and  Structure  Functions 

The  follow  SDF  form  is  often  used  because  it  admits  both  a  power-law  regime  and  an 
analytic  inverse[41] 


Q(oj)  =  C 


r>/2) 


\A  +  v/2k> 


7V1 

Qs 


(5.40) 


We  will  now  use  this  form  to  develop  a  relation  between  the  structure  constant  and  the 
spectral  strength  parameter.  For  small  arguments  Kl/(x)  &  (1/2)  F(z/)  |x/2|~^L  Thus,  as 
long  as  qi  «  u  <<  qs  and  qi/qs  <<  1  (5.40)  has  the  desired  power-law  form 

Q(cj)  x  Cu~v .  (5.41) 

This  “self-similar”  Fourier  transform  is 

R(v)  =  ^2(^-1)/2f(^72)  P2  +  +  (QlvW/qT1’  (5-42) 


where  R(y)  =  E  [x(t)x*  (t1)]  and 


r  =  qilqs- 

Unlike  the  SDF,  the  covariance  function  does  admit  a  scale  free  power-law  form, 
that, 


E[\x\2)  =  R( 0)  » 


c  r((»/- 1)/2) 

>/5F  T(u/2 )qrl  ’ 


Consider 
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which  is  dominated  by  the  outer-scale  wavenumber  through  the  factor  qlL~v . 

The  coherence  properties  of  power-law  process  with  v  >  1  are  better  characterized  by 
the  structure  function 


D(y) 


x 


=  E[\x(t)-x(f)\*]  =  2{R{0)-R{y)) 

_  2C  iV 

“  v^F  2Cv-i>/2r(«//2) 


\Jr~  +  (^Lj/)2'1'  1)/2/£>-i)/2(\A2  +  ( qtW)/N 


K 


V-  1 


(5.43) 

(5.44) 


where 

AT  =  r(*'+i)  Ar(t/_1)/2(r)  «  r((i/  -  l)/2)2^-1^2.  (5.45) 

The  term  in  square  brackets  converges  to  zero  as  qi  approaches  zero.  Thus,  an  application 
of  LHopital’s  rule  yields  the  asymptotic  approximation 


D(y) 


C  r((3— v)/2)23~ 
nA  (u-i)r(t//2> 

2  l  t*'— 1 

=  ^  M 


■M 


i/— i 


(5.46) 
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SpecLab-Scripts 


A-lFilt_Par 


Prediction  of  a  process  (s) . 

Model: 

y=s+n=sm+sr+n 
n  -  white  noise  (not  stationary) 
s=sm+sr  with  sm  (constant)  mean  value 

Usage: 


function 

SP: 

=F ilt _Par ( ind , y , vs , vn , 1) 

Input : 

y 

- 

the  observations 

column  vector 

ind 

argument  vector  (in  sampling  units) 

n 

vs 

variance  of  sr  (sr  being  zero  mean) 

scalar 

vn 

- 

variance  of  ’measurement  noise’  n 

column  vector 

1 

correlation  range  of  sr  in  sampling  units 

scalar 

(sr 

stationary,  n  white  not  stationary) 

Output 

:  sp 

- 

filtered  version  of  y  (prediction  of  s) 

A-2GetParproc 


Estimates  and  returns  the  power  law  parameter  processes 

given  the  segmentation,  the  scale  spectra  in  the  (spatial)  segments 

and  the  estimates  of  the  inertial  range. 

Usage 

P=ScaleSpec(Z) 

Input 

Z  -  Matlab  data  structure  containing  spectra 

and  inertial  ranges.  Initialized  by  InitSpec, 

ScaleSpec  and  InertialSpec  (see  also  InitZ) . 

Output 

P  -  The  estimated  parameter  process 

See  Also:  ScaleSpec  ScfiltSpec  SynthSpec 


A-3GetVario 


The  variogram  of  the  input  process  is  estimated  interactively. 
Usage 

[intcept  1  v] =GetVario (y, titles) 
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Input 

y  -  the  input  data  vector 

titles  -  title  string 

Output 

1  -  estimate  of  the  correlation  range 

v  -  estimate  of  the  variance 

intcept  -  estimate  of  the  intercept 

See  Also:  VarioSpec  VarioEst 


A~3Get  Jnertial 


Estimates  the  inertial  ranges  given  the  scale  spectra 
in  the  set  of  segments  and  the  length  of  the  data  vector 
from  which  the  scale  spectra  have  been  computed. 

The  estimation  is  done  by  comparing  deviation  of  the  log  scale 
spectrum  from  a  linear  log  spectrum  with  the  mean  value 
of  this  measure  for  Brownian  motion. 

Usage 

[inertial]  =Get_Inertial (scspect , dyad) ; 

Inputs 

scspect  -  the  set  of  scale  spectra  (#seg  by  #scales  vector) 
dyad  -  dyadic  dimension  of  original  data  vector 

Outputs 

inertial  -  the  inertial  ranges  (#seg  by  2  vector) 

See  Also:  InertialSpec 

This  matrix  is  stored  on  file:  Ref res 
and  contains  the  measure  of  deviation 
for  Brownian  motion. 

Dimensions  is  maxdyad  by  maxdyad  with  first  index 
corresponding  to  dyadic  dimension  and  second  index  to 
size  of  inertial  range. 

Minimum  number  of  scales  in  inertial  range  (>  1) 

A-3HelpSpec 


Ref  res  7. 


mins  7, 


Gives  online  help  for  the  main  menu  SpecLab  options 
Usage 

HelpSpec 
See  Also:  Main 
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The  fields  of  Z  are: 

data  -  the  data  vector 

wavec  &  wavel  -  c  &  1  from  wavedec  (length  of  data  is  1 (length (1)  ) 
nseg  -  number  of  (uniform)  segments 

in  segmentation  (scalar) 
scspect  -  scale  spectrum 

(nseg  by  m  array;  [short  to  long  scales]; 
note  that  m  is  number  of  scales  stored  in 
wavec,  which  is  ‘ length (wavel) -2 7 ) . 
inertial  -  inner  and  outer  scales  in  the  segments 

(nseg  by  2  array;  [inner  outer]) 

parproc  -  the  exponent  and  intercept  estimates  in  segments 

(nseg  by  2  array;  [H  interc] ) 
filtpar  -  parameters  used  in  filtering  of  parproc 

(1  by  6  array; 

[intceptH  varH  corrlengthH  . . . 
intceptLoglnt  varLoglnt  corrlengthLoglnt] ) 
filtparproc  -  filtered  parameter  process 

trueparproc  -  true  parameter  process  (if  process  is  simulated) 

See  Also:  SpecLab  Main 


A-3Linreg 

Estimates  the  power  law  parameters  for  the  scale  spectrum. 
Usage 

[inspect  ,c,H]  =Linreg(spect , low, high) 

Input 

spect 
low  high 

Output 

mspect 
c 
H 

See  Also:  ProcSpec 


A- 3  Main 

Main  is  part  of  SpecLab’ s  GUI  driven  interface  that  allows  the  user 
to  apply  the  components  of  SpecLab. 

Main  presents  the  user  with  the  main  menu  that  allow  him  or  her 
to  apply  the  different  estimation  or  plotting  capabilities 
of  SpecLab.  These  are  implemented  by  intermediate  level 
routines  that  call  SpecLab 7 s  basic  or  tool  routines. 


-  scale  spectrum  (column  vector) 

-  inertial  range  used  in  estimation 

-  estimated  model  spectrum 

-  estimated  log  intercept 

-  estimated  Hurst  exponent 
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Z=Main(Z) 


Usage 


Output 

Z  -  Matlab  data  structure  conaining  the  analyzed 
process  and  its  spectral  information 
(see  InitZ  for  format) . 

See  Also:  SpecLab  Z  HelpSpec  InertialSpec  ProcSpec  ScaleSpec 

ScfiltSpec  SmoothSpec  SynthSpec  VarioSpec 


A“3MakeRefres 


Computes  the  percentail  of  the  residual  of  the  scale  spectrum  relative  to 
the  estimated  linear  model  for  various  dyadic  lengths  of  data 
and  sizes  of  the  scale  range.  Result  is  used  when  inertial  range  is 
estimated . 

Usage 

MakeRef res 


Outputs 

The  matrix  Refres  is  stored  on  the  file  Refres 

with  storage  format  (ascii) .  Refres  is  a  maxdyad  by  maxdyad 

array  reshaped  to  linear  form  for  storage. 

See  Also:  Get_Inertial 

Parameter 


A-3ProcSpec 


Estimates  plots  and  returns  the  power  law  parameter  processes 
given  the  segmentation,  the  scale  spectra  in  the  (spatial)  segments 
and  the  estimates  of  the  inertial  range. 


Usage 

Z=ScaleSpec (Z) 

Input 

Z  -  Matlab  data  structure  containing  spectra 

and  inertial  ranges.  Initialized  by  InitSpec, 
ScaleSpec  and  InertialSpec  (see  also  InitZ) . 


Output 

Z  -  The  input  Matlab  data  structure  with 

the  parameter  processes  added. 


See  Also :  Main 
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A-3ReadSpec 


Reads  data  vector 
Usage 

[Z ,  path.]  =ReadSpec  (mode) 

Inputs 

mode  -  Type  of  data  input 

1:  binary 
2:  ascii 


Output 

Z  -  Speclab  data  structure 

path  -  path  is  0  if  no  file  was  read 


See  Also:  SpecLab  initZ  getdat 


A-3ScaleSpec 


Computes  plots  and  returns  the  wavelet  scale  spectrum 
given  the  data  and  its  wavelet  decomposition.  The 
user  can  choose  the  number  of  segments  in  a  uniform 
segmentation  of  the  data. 


Usage 

Z=ScaleSpec (Z) 

Input 

Z  -  Matlab  data  structure  containing  data  vector  and 

its  wavelet  decomposition.  Initialized 
by  InitSpec  (see  also  InitZ)  . 

Output 

Z  -  Matlab  data  structure  that  contains  scale  filtered 

data,  its  wavelet  transform,  the  chosen  segmentation 
and  scale  spectra  in  segments. 

See  Also  Main 


A-3Scale_Filt 


Low  pass  scale  filtering  of  data 
Usage 

[X,c,l]=scale_filt(c,l, inner) ; 
Inputs 


62 


c,l  -  wavelet  transform  obtained  from  wavedec 

inner  -  the  inner  scale  used  in  filtering 

Outputs 

c,l  -  low  pass  truncated  wavelet  transform 

X  -  low  pass  filtered  data 

See  Also:  ScfiltSpec 


A-3ScfiltSpec 


Scale  filters  data,  plots  them 

and  returns  data-structure  with  spectral  information. 

Usage 

Z=InitSpec (Z) 

Inputs 

Z  -  Data  structure  containing  spectral  information 

(see  InitZ  for  format) . 

Output 

Z  -  Matlab  data  structure  that  contains  a  scale  filtered 

version  of  input  data. 

See  Also:  Main 


A-3SegSim 


Simulates  a  locally  stationary  power  law 
Usage 

[z] =SegSim(dyad ,pv , iv) 

Inputs 

dyad  -  dyadic  dimension  of  simulation  grid 

pv  -  Hurst  process 

iv  -  log  intercept  process 


Outputs 

z  -  the  simulated  process  (2~dyad  by  1  vector) 
See  Also:  SynthSpec 


A-3SimExp 
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Simulates  the  locally  stationary  power  law  parameters. 

The  model  is  a  truncated  Gaussian  exponentially  correlated 
stochastic  process. 

Usage 

[pv, iv] =SimExp(dyad,dyadl,Hm,Hs ,1s ,mode) 

Inputs 

dyad  -  dyadic  dimension  of  parameter  process 
dyadl  -  dyadic  dimension  of  corrlation  length 
Hm  -  mean  of  Hurst  parameter 

Hs  -  relative  fluctuations  in  Hurst  parameter 

Is  -  relative  fluctuations  in  log  intercept 

mode  -  optinal  parameter 

mode  :  1  -  Triangular 

2  -  Gaussian  correlation 

3  -  Exponential  correlation 

Outputs 

pv  -  Hurst  process 

iv  -  log  intercept  process 

See  Also:  SynthSpec 


A-3SmoothSpec 

The  power  law  parameter  estimates  are  smoothed  according 
to  the  model  estimated  for  their  stochastic  variations. 

The  smoothed  are  plotted  optionally  with  the  unsmoothed 
versions . 

Usage 

Z=SmoothSpec (Z) 

Input 

Z  -  Matlab  data  structure  containing 

the  (unfiltered)  parameter  processes  and 
variogram  parameters  for  their  spatial  variation, 
(see  also  InitZ) . 

Output 

Z  -  The  input  Matlab  data  structure  with 

the  smoothed  parameter  processes  added. 

See  Also:  Main 


A-3SpecLab 

SpecLab  is  an  GUI  driven  interface  that  allows  the  user 
to  apply  the  components  of  SpecLab  to  a  set  of  test 
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data  sets  or  to  simulated  data  sets  whose  parameters  are 
determined  by  the  user.  The  user  can  alternatively  load 
his  own  dataset. 

Using  the  GUI  interface  the  user  can  estimate  a  Local  Power 
Law  model  for  the  data  and  also  plot  the  dataset  or 
its  spectrum  in  various  ways. 

Usage 

Z=SpecLab 

Output 

Z  -  Matlab  data  structure  conaining  the  analyzed 
process  and  its  spectral  information 
(see  InitZ  for  format) . 

See  Also  InitZ  Main 


A-3Spect 


Computes  the  wavelet  scale  spectrum 

given  the  wavelet  decomposition  and  the  number  of 

segments  in  a  uniform  segmentation  of  the  data. 

Usage 

[scspect] =Spect (c ,1 ,nseg) 

Input 

c  &  1  -  c  &  1  from  wavedec 

nseg  -  number  of  (uniform)  segments 

in  segmentation  (scalar) 

Output 

scspect  -  scale  spectrum 

(nseg  by  m  array;  [short  to  long  scales] ; 
note  that  m  is  number  of  scales  stored  in 
c,  which  is  ‘ length (1) -2 } ) . 

See  Also:  ScaleSpec  ScfiltSpec  SynthSpec 


A-3SynthSpec 


Simulates  a  locally  stationary  power  law. 

The  user  specify  the  power  law  parameters  from  a  simple  model. 
The  intercept  and  Hurst  exponent  of  the  power  law  are  taken  to 
be  realizations  of  a  stochastic  process. 

Also  inner  and  outer  scales  are  specified. 

First  a  pure  power  law  is  simulated.  Information 
in  scales  beyond  the  inertial  range  is  then  tapered. 
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Usage 

[Z] =SynthSpec 
Output 

Z  -  Matlab  data  structure  containing  the  simulated  data 

vector  and  its  spectral  information 
(see  InitZ  for  information  about  the  fields). 

See  Also:  main 


A-3Taper 

Tapers  the  locally  stationary  power  law  according  to 
inner  and  outer  scales  with  Gaussian  tapering. 

Usage 

[z] =Taper (z , ins , outs) 

Inputs 

z  -  the  process 
ins  -  the  inner  scale 
outs  -  the  outer  scale 


Outputs 

z  -  the  tapered  process 
See  Also:  SynthSpec 

A-3VarioEst 


Computes  the  variogram  of  the  input  process 

and  fit  it  to  a  vertically  shifted  exponential  model  using  least  squares. 
Usage 

[intcept  1  v  vario] =VarioEst (y ,maxlag) ; 

Inputs 

y  -  process  for  which  the  variogram  is  computed 

maxlag  -  the  maximum  lag  for  which  the  variogram  is  computed 
Outputs 

intcept  -  estimated  intercept 
1  -  estimated  correlation  length 

v  "  variance 

The  fitted  variogram  model  is:  intcept+v* (l-exp(-x/l) ) 

vario  -  the  empirical  variogram  of  y 


See  Also:  GetVario  regmodel 


A-3VarioSpec 
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Estimates  the  variogram  or  structure  function  of  the 
the  power  law  parameter  processes.  The  estimate 
can  interactively  be  modified  by  the  user. 

Usage 

[Z] =VarioSpec (Z) ; 

Input 

Z  -  Matlab  data  structure  containing 

the  (unfiltered)  parameter  processes, 
(see  also  InitZ) . 

Output 

Z  -  The  input  Matlab  data  structure  with 

the  variogram  parameters  added. 

See  Also  Main 


A-3WriteSpec 

Reads  data  vector 
Usage 

WriteSpec (Z) 

Inputs 

Z  -  Matlab  data  structure  with  process 
(see  InitZ  for  format) 

See  Also:  main  InitZ  writedat 


A- 3 aero 


This  workout  illustrates  how  the  figures  in  the  paper: 

Wavelet  based  estimation  of  local  Kolmogorov  turbulence 
by  George  Papanicolaou  and 
Knut  Solna 

can  be  generated  using  SpecLab.  It  also  serves  as  an  introduction 
to  SpecLab  itself. 

Usage 

Z=aero 

Output :  Z  -  data  structure  with  aerothermal  data 

(see  InitZ  for  format) 

See  Also:  SpecLab  InitZ 
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A-3displayP 


Plots  the  power  law  parameter  processes 
Usage 

displayP(Z,mode) 

Inputs 

Z  -  Data  structure  containing  data  and  spectral  information 

mode  -  type  of  plot 

1  :  Plot  unsmoothed  parameters 

2  :  Plot  smoothed  and  unsmoothed  parameters 

3  :  Plot  smoothed  parameters 

See  Also:  Main  InitZ  ProcSpec  SmoothSpec 


A-3display3 

Plots  the  set  of  scale  spectra 
Usage 

displays (scspect ,N, mode, inertial) 

Inputs 

scspect  -  scale  spectrum 

(nseg  by  m  array;  [short  to  long  scales]; 
with  m  the  number  of  scales  and  nseg  the 
number  of  spatial  segments. 

N  -  length  of  data  vector  on  which  the  scale 

spectrum  is  based, 
mode  -  type  of  plot : 

1  -  color  image 

2  -  graphs 

3  -  color  image  with  inertial  range 

4  -  graphs  with  inertial  range 

inertial  -  inner  and  outer  scales  in  the  segments 

(nseg  by  2  array;  [inner  outer]) 

See  Also:  Main  InertialSpec  ScaleSpec  SynthSpec 


A-3displaySyn 

Plots  simulated  process  and  parameters 
Usage 

displaySyn(X,pv, iv) 

Inputs 

X  -  simulated  process 

pv  -  the  Hurst  process 
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IV 


-  the  log  intercept  process 


See  Also:  Main  SynthSpec 


A-3displayVario 


Plots  variogram 
Usage 

display Var io (model , var io , titles) 

Inputs 

model  -  model  variogram 

vario  -  empirical  variogram 

intcept  -  the  itercept  of  the  variogram 

titles  -  the  title  string 

See  Also:  Main  GetVario 


A-3displayW 


Plots  level  of  detail  coefficients 
Usage 

displayW(c , 1) 

Inputs 

c,l  -  Wavelet  decomposition 


See  Also:  Main  ScfiltSpec 


A-3displayX 


Plots  data  vector 
Usage 

displayX(X , titlestr) 

Inputs 

X  -  The  data 

titlestr  -  title  of  plot 

See  Also:  Main  ScfiltSpec 


A-3getdat 


Reads  data  vector 
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Usage 

[Z,path]=getdat (f ilen, path, mind) 

Inputs 

filen  -  file  name 

path  -  path 

mind  -  minimum  dyadic  length  of  data-vector 

mode 

1:  binary 

2:  ascii 

Output 

Z  -  Data  vector 

path  -  path  is  0  if  no  vector  was  read 

See  Also  ReadSpec  InitZ 


A-3hfunc 

Computes  the  function  h(H)  that  is  involved  in  the  definition  of 
the  mean  and  variance  of  scale  spectral  points. 

Usage 

h=hfunc(H) 

Input :  H  -  the  Hurst  exponent 

Output  h(H)  -  scaling  factor  for  mean  and  variance 

of  scalespectral  points. 

See  Also:  Linreg 


A~3regmodel 

Gives  model  variogram  or  structure  function 
Usage 

y=regmodel (beta,x) ; 

Inputs 

beta  -  correlation  length 

x  ~  lag 

Outputs 

y  -  value  of  variogram 

See  Also:  VarioEst 


A-3showinput 
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Shows  the  input  menu  with  a  paricular  entry  marked 
Usage 

showinput (opt) 

Input : 

opt  -  the  menu  entry  that  is  to  be  illustrated  next 
See  Also:  aero 

A-3sho\vmenu 

Shows  the  main  menu  with  a  paricular  entry  marked 

Usage 

showmenu(opt) 

Input : 

opt  -  the  menu  entry  that  is  to  be  illustrated  next 
See  Also:  aero 

A-3sim_ob 


Simulates  OBSERVED  fBm 

The  process  is  scaled  so  that  the  value  of  the 
structure  function  at  maximum  lag  (lag=2~d)  is  1 
and  condition  on  the  process  being  zero  at 
the  ’ zero' th  value  of  the  index. 

The  measurement  model  corresponds  to  the  observations 
being  integrals  of  the  fBm  over  unit  intervals. 

Input  :  d  -  dyadic  dimension 

H  -  Hurst  exponent ! 

Output:  xold  -  Simulated  sequence  (2~d  by  1  vector) 

Parameters:  D  -  dyadic  dimension  of  supergrid  in  simulation 
1  -  parameter  that  determines  the  SIZE 

of  the  conditioning  set  in  the  simulation, 
(simulation  ACCURACY  and  also  TIME 
increases  with  the  value  of  these 
two  parameters) .  Note  paremeters 
need  to  satisfy  2~D  >  2*1  >  2. 
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A-3simcoef 


Computes  coefficients  used  by  sim_ob  in  the  simulation 


of  observed  fBm. 

Input  : 

nl  ,n2 
del 

H 

-  dimensions  of  conditioning  sets 

-  lag  between  conditioning  nodes 

-  the  Hurst  coefficient 

Output : 

cl ,  c2 

s 

-  cl,c2  coefficients  of  conditioning  nodes 

-  conditional  standard  deviation  of  simulation 

node . 

See  Also 

:  sim_ob 

A-3vfunc 


Computes  the  value  of  the  structure 
function  for  fBm  assuming  a  measurement 
model  where  integrals  of  B  over  unit  intervals 
are  observed,  (value  computed  assuming  \sigma’'2=l) 

Useage:  function  v=vfunc (lag,H) 


Input : 


Output 


H  -  the  Hurst  exponent 

lag  -  separation  distance  at  which  structure 

function  is  computed  (NB  an  integer) 
v  -  value  of  structure  function 


See  Also:  simcoef  sim_ob 


A-3writedat 


Writes  data  vector 
Usage 

[path] =wr itedat (f ilen , path , Z) 

Inputs 

filen  -  file  name 

path  -  path 

Z  -  Matlab  data  structure  with  process 

(see  InitZ  for  format) 


Output 

path  -  path=0  if  file  not  propely  written 
See  Also:  WriteSpec  InitZ 
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